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19.ABSTRACT 

^^>-This  report  describes  tho-restths  of  atwcryear  research  projectgafr  continuum  modeling 
and  vibration  control  of  flexible  structures  with  application  to  active  control  of  vibrations 
in  large  space  structures.  A  comprehensive  methodology  is  discussed  for  the  construction 
of  effective  (linear')  modelg  for  lew^e  composite  structures  consisting  of  various  flexible 
members  (e.g.  beams,  trusses,  etc.)  and  rigid  body  elements.:  Since  our  ultimate  concern 
is  the  active,  feedback  control  of  such  systems,  we  fin  dwt 'convenient  to  concentrate  on 
frequen^j^dgmain  modeling. -yWe  start  at  the  component  level  and  shp^a  systematic 
procedureAlor  computing  the  irrational  transfer  functions.  #equiiClf. '  Then  by  standard 
transform  methods  a  complete  hybrid  model  is  developed.  The  methods  were  coded  in  a 
computer  algebra  system  (SMP  was  used)  which  automated  the  model  building  process 
and  produced  Fortran  code  for  numerical  evaluation  of  the  frequency  responses.  / 

Under  certain  conditions  the  dynamics  of  such  large,  low  mass  structures^  having  a 
regular  infrastructure,  can  be  modeled  by  the  dynamics  of  continuua;  e.g.,  trusses  can 
be  modeled  as  elastic  beams.  We  demonstrate  hov^ffective  continuum  models  of  lattice 
structures  with  regular  infrastructure  can  be  obtained  by  a  systematic  procedure  based  on 
an  asymptotic  analysis  of  multiple  scales  called  homogenizationf^This  method  is  applied 
to  several  examples  and  it  is  jghown  tha{  accurate  computation^?  the  required  parameters 
of  such  continuum  models  ts^somewhat  more  subtle  than  merely  averaging  over  lattice 
cells.  \Certain  “corrector”  formulae,  required  for  such  computations,  are  derived  for  several 
examples.  In  addition  the  nature  of  the  process  of  combining  homogenization  and  optimal 
control  is  discussed  in  detail  _ 

For  the  computation  of  distributed  parameter  control'we  concentratetTdn^fcn  optimal 
frequency  domain  method^  based  on  solving  an  associated  Wiener ^Hopf  problem.  The 
method'which  was  pOptrlarized'by  ^OirDavisl*toploys  effective  numerical  algorithms  (e.g. 
FFT,  etc.)  to  compute  a  certain  spectral  factorization  of  a  possibly  matrixrvalued  (in  the 
multiple  control  case)  Hermittian,  positive-definite  transform  by  sampling  the  frequency 
response.  The  control  laws  considered  in  this  report  take  the  form  of  distributed  state 
feedback  with  respect  to  a  naturally  defined,  distributed  statejkpace  of  functions  over  the 
spatial  domain  of  the  structure.  - 


Several  examples  are  considered  in  detail  and  certain  numerical  sensitivities  are  noted 
for  the  method.  The  relationship  between  this  method  for  distributed  control  and  more 
standard  methods  based  on  finite-dimensional,  reduced-order  modeling  is  discussed.  In  this 
context  the  distributed  parameter  method  raises  several  basic  questions  with  respect  to 
continuum  modeling  which  never  arise  when  finite-dimensional  approximations  are  used 
from  the  outset.  In  particular,  as  discussed  by  Gibson,  the  nature  of  various  standard 
models  for  structural  damping  is  deemed  crucial. 

Our  ultimate  goal  in  this  project  was  to  demonstrate  both  effective  methods  for  con¬ 
tinuum  modeling  of  low  mass  density  structures  and  computation  of  stabilizing,  optimal 
control  laws  for  these  systems.  Examples  indicate  that  stabilization  and  control  of  essen¬ 
tially  all  controllable  modes  of  the  given  system  will  be  difficult  to  achieve  in  practice  using 
a  finite  array  of  sensors.  However,  it  is  shown  that,  under  certain  conditions,  an  aribitrary, 
finite  number  of  modes  can  be  controlled  using  these  methods  without  destabilizing  the 
remaining  uncontrolled  modes. 
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with  Application  to  Vibration  Control” 

to  be  submitted  to  AIAA  Journal,  and 

W.H.  Bennett,  H.G.  Kwatny,  and  T.S.  Aven,  “Practical  Issues  in  Computation 
of  Optimal  Distributed  State  Feedback  Control  for  Flexible  Structures” 

to  be  submitted  to  AIAA  Guidance,  Navigation,  and  Control  Conference,  August  1987, 
Monterey,  California. 
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Executive  Summary 


In  this  project  we  considered  modeling  the  flexible  dynamics  of  large  space  structures 
for  the  purpose  of  designing  active  control  laws  for  the  suppression  of  vibrations.  We 
focus  on  distributed  parameter  models  which  capture  the  spatial  dependency  of  the  in¬ 
teraction  between  control  forces,  load  forces  (on  board  systems),  exogenous  distrubances, 
and  sensors.  Standard  methods  for  such  analysis  used  throughout  the  aerospace  industry 
are  based  on  finite-element  analysis.  In  contrast  in  this  project  we  focused  on  the  use  of 
continuum  models  for  the  elastic  response  of  certain  mechanical  components  like  beams, 
cables,  membranes.  A  composite  model  for  a  space  structure  may  consist  of  several  of 
these  elastic  elements  interacting  at  their  boundaries  with  rigid  body  elements.  Effective 
continuum  models  for  periodic  truss  structures  may  also  be  used  to  simplify  the  model 
construction.  We  have  employed  an  asymptotic  method  of  mutliple  time  and  spatial  scales 
called  homogenization  to  compute  such  continuum  models  for  a  variety  of  truss  structures 
under  several  assumptions. 

Control  for  distributed  systems  centers  on  the  definition  of  a  distributed  state  for  the 
mechanical  system.  A  method  for  computation  of  a  distributed  feedback  control  gain  is 
developed  in  detail  including  numerical  algorithms  and  considerable  testing  is  included. 

In  the  first  section  we  provide  a  technical  introduction  to  the  problem  of  control  of 
distributed  parameter  systems  arising  from  elastic  mechanical  structures.  The  second 
section  contains  a  detailed  discussion  of  a  relatively  new  approach  to  computing  distributed 
control  gains  for  such  systems.  The  method  is  based  on  a  classical  Wiener-Hopf  problem 
and  involves  considerable  numerical  computations.  We  provide  comparison  with  standard 
modal  and  finite-element  methods  and  consider  the  practical  limitation  of  the  approach. 

The  third  section  contains  a  discussion  of  frequency  response  computations  for  con¬ 
tinuum  models  arising  from  elastic  dynamics  modeled  by  partial  differential  equations. 
Composite  structures  are  also  considered  in  this  section  and  it  is  shown  that  standard 
computations  based  on  Laplace  transforms  can  be  used  effectively  to  computed  the  re¬ 
quired  transform  for  rather  complex  interconnections  between  elastic  and  rigid  compo¬ 
nents.  In  section  four  we  consider  homogenization  of  regular  structures.  The  method 
is  demonstrated  on  a  simple  one  dimensional  example  which  shows  the  existence  of  cer¬ 
tain  ‘corrector’  formulae  which  are  important  in  computation  of  such  models.  A  major 
question  is  the  effectiveness  of  control  laws  developed  for  the  resulting  continuum  model 
when  applied  to  the  finite  structural  lattice.  In  this  section  this  question  of  combined 
homogenization  and  control  is  considered  in  some  detail. 

Finally  section  three  contains  several  examples  which  illustrate  the  combined  modeling 
and  control  methods  proposed. 

The  problem  of  deriving  complete  transfer  function  representations  of  the  structural 
models  including  transient  response  characterization  by  Green’  functions  is  completely 
solved  in  this  work.  We  show  that  the  analysis  of  certain  types  structures — beams  with 
one  or  more  degrees  of  freedom — can  be  “automated”  (though  this  is  far  from  trivial) 
using  a  computer-based,  symbolic  manipulation  system.  The  major  question  remains  the 
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numerical  evaluation  of  these  functions.  It  is  possible  to  use  a  variety  of  expansions  and 
alternate  algebraic  representations  for  the  functions  for  various  ranges  of  their  domains. 
This  analysis  can  be  supported  with  symbolic  computation.  Also,  it  appears  likely  that 
a  truely  functional  methodology  will  require  capability  for  rational  approximation  and 
interpolation  of  these  functions.  These  issues  limit  the  practical  application  of  the  method 
in  question.  Computer  code  for  the  symbolic  and  numerical  algorithms  are  contained  in 
the  appendices. 
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Introduction 


It  is  now  generally  accepted  that  large,  low  mass  density  lattice  structures  will  be  essential 
for  several  near  term  space  applications.  Moreover,  it  is  apparent  that  active  control  of 
structural  vibrations  will  be  necessary  to  enhance  their  stiffness  and  damping  properties. 
In  this  report  we  consider  the  construction  of  effective  mathematical  models  for  elastic 
dynamics  of  space  structures  with  the  objective  of  designing  active  control  laws  for  these 
systems. 

The  success  of  active  control  for  such  structures  will  hindge  to  some  extent  on  the 
ability  of  a  control  law  to  react  to  vibratory  responses  which  may  be  initially  localized 
before  they  propagate  throughout  a  structure.  This  leads  naturally  to  questions  of  how 
to  implement  active  control  so  as  to  distribute  the  control  effort  spatially  as  it  is  needed. 
We  argue  that  well  known  methods  exist  for  the  control  of  distributed  parameter  systems 
and  can  be  effectively  applied  if  continuum  models  for  the  candidate  space  structures  can 
be  computed.  The  nature  of  the  required  models  is  however  quite  different  from  the  more 
standard  finite  element  models  which  are  popular  for  large  structural  analysis  problems 
throughout  the  aerospace  industry. 

We  begin  in  this  section  with  a  review  of  continuum  models  for  active  structural  con¬ 
trol.  We  highlight  the  nature  of  abstract  state  space  models  for  these  systems.  In  this 
study  we  have  employed  one  method  for  the  computation  of  a  distributed  control  law 
for  continuum  dynamics  based  on  a  numerical  procedure  for  spectral  factoriztion.  This 
method  requires  the  computation  of  certain  representation  for  an  underlying  state-space 
model  for  the  structure  to  be  controlled.  In  a  later  section  we  discuss  the  theoretical  basis 
for  computing  effective  distributed  parameter  models  for  large  truss  structures  with  ran¬ 
dom  lattice  infastructure.  The  method  which  involves  “homogenization”  (an  asymptotic 
analysis  of  multiple  scales)  leads  to  the  well  known  Timeshenko  model  for  beam  dynamics. 
The  analysis  provides  formulae  for  the  effective  beam  parameters  which  are  quite  different 
than  have  been  suggested  by  other  averaging  schemes  [1,2,3]. 

Comprehensive  models  of  flexible  spacecraft  dynamics  will  involve  systems  with  fairly 
complex  interconnections  of  lumped  and  distributed  subsystems,  and  therefore,  we  intend 
to  construct  the  overall  models  by  first  developing  subsystem  models  and  then  combining 
them  according  to  the  required  interconnection  rules.  These  interconnections  lead  to  basic 
questions  of  causality  and  well-posedness  of  certain  standard  models  for  beams.  These 
questions  are  crucial  to  the  computation  of  hybrid,  state-space  modeling  of  an  integrated 
space  platform. 

Throughout  this  effort  we  have  focused  on  the  potential  for  automatic  and  computer- 
aided  computation  of  the  models  by  a  combination  of  modern  computer  algebra  [4]  (sym¬ 
bolic  manipulation)  and  numerical  methods.  In  our  efforts  we  have  used  the  program  SMP 
[5,6].  We  will  review  the  progress  in  using  SMP  for  the  computation  of  irrational  transfer 
function  models  for  hybrid  problems  in  a  later  section. 
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1.1  Generic  Models  for  Structural  Dynamics 

In  this  section,  we  consider  a  generic  model  for  elastic  dynamics  of  structures  from  the 
point  of  view  of  continuum  modeling.  We  will  summarize  the  construction  of  a  state  space 
model  and  introduce  a  typical  control  problem  for  vibration  suppression.  We  highlight  the 
modal  approximations  which  are  popular  for  these  problems  and  proceed  to  demonstrate 
an  effective  alternate  technique  for  model  construction  and  control  computation  based 
on  the  semi-group  property  [8]  of  a  state  space  model.  Effectively,  modeling  and  control 
law  computation  can  proceed  in  the  frequency  domain,  based  on  transfer  function  meth¬ 
ods,  permitting  the  direct  computation  of  a  resolvent  operator.  We  focus  on  the  class  of 
structural  control  problems  for  which  the  question  of  control  of  propagation  of  wave-like 
disturbances  is  important.  In  this  framework,  we  can  present  the  semi-group  theory  by 
concrete  computations  of  practical  interest  to  structural  and  control  system  engineers. 

A  popular  continuum  model  for  a  flexible  structure  [9]  is  described  by  a  system  of 
partial  differential  equations  (PDE) 


m{z) 


d2w{t,z)  ,  „  dw(t,z) 

- 57; —  +  — 


+  A0w(t,z)  =  F[t,z) 


where  tv(t,z)  is  an  N-vector  of  displacements  of  a  structure  f 2  with  respect  to  some  equi¬ 
librium  for  H  is  a  bounded,  open  set  in  l.  The  (vector)  z  6  0  is  a  coordinate  in  fi. 
We  assume  the  boundary  dfl  is  smooth.  The  mass  density  m(z)  is  positive  definite  and 
bounded  on  dfl.  The  damping  term  Dodw/dt  models  both  (asymmetric)  gyroscopic  and 
(symmetric)  structural  damping  effects.  The  internal  restoring  force  A0 w  is  generated  by 
a  time-invariant,  differential  operator  A0  for  the  structure.  For  most  common  structural 
models,  A0  is  an  unbounded,  differential  operator  with  domain  D(A0)  consisting  of  cer¬ 
tain  smooth  functions  satisfying  appropriate  boundary  conditions  on  dfl.  Thus,  for  these 
problems,  D(A0)  is  typically  dense  in  the  Hilbert  space  =  L2(fl)  endowed  with  its  nat¬ 
ural  inner  product  (x,y)0  =  /n  xr(2)y(z)  dz.  Often  (but  not  always),  the  spectrum  of  Ao, 
<r(A0),  consists  of  discrete  eigenvalues  with  associated  eigenfunctions  which  constitute  a 
basis  for  £2(n)2  . 

The  applied  force  distribution  F(f,z)  can  be  thought  of  as  consisting  of  three  compo¬ 
nents 

F{t,z)  =  Frf(f,z)  +  Fe{t,z)  +  Fa{t,z)  (2) 

where  Fi  is  N-vector  of  exogenous  disturbances  (possibly  forces  and  torques),  Fe  is  a 
continuous,  distributed  controlled  force  field  (an  available  option  in  only  some  special 
applications),  and  F„  represents  controlled  forces  due  to  localized  actuation; 

F*{t,z)  =  £M«)*,(f)  =  -Bou(f).  (3) 

>=i 

The  actuator  influence  functions  6;(z)  are  highly  localized  in  f)  and  can  be  approximated 
by  Dirac  delta  functions.  We  assume  that  a  finite  number  p  of  measurements  can  be  made 

*R  denote*  the  reel  number*. 

aT*  ia  the  apace  of  function*  *uch  that  /  €  £p(0)  then  (/n  |/(w)|pdu;),/a 
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y{t)  =  C0w  +  C'0^j-  (4) 

where  y(t)  is  a  p- vector.  The  operators  B0  :  SRm  — ►  Ho,  Co  ■  Ho  — *  3RP,  and  C'0  :  Ho  — »  5RP 
are  bounded. 

The  standard  vibration  control  problem  for  this  model  is  to  find  the  controls  u;(f),  j  = 
1, . . . ,  k  (we  ignore  the  possibility  of  Fe)  given  the  observations  y(t)  to  maintain  the  system 
state,  e.g., 


,  (w(t,  *)\ 

’ J  Um)/ 


as  close  to  its  equilibrium  state  as  possible. 


1.1.1  State  Space  Models 

The  choice  of  state  space  given  by  (5)  is  often  attractive  for  models  in  the  generic  form 

(l).  (We  will  discuss  later  that  attractive  alternate  state  space  models  can  arise  in  hybrid 

constructions.)  A  natural  assumption  for  structural  problems  is  that  Aq  is  symmetric  with 

compact  resolvent  and  discrete  (real)  spectrum  which  is  bounded  from  below  [9].  The 

1 

state  (5)  can  then  be  considered  as  an  element  of  a  Hilbert  space  H  =  D(Ajj)  x  Ho  with 
the  energy  norm 

||z|||  =  {«>,  Aow)0  +  (mw,  w)0  (6) 

where  the  first  term  represents  potential  energy  and  the  second  term  is  kinetic  energy. 
Thus  an  (abstract)  state  space  model  can  be  written 

i(t,z)  =  Ax(t,z)  +  Bu(t)  (7) 

y(t)  =  Cx(t,z) 

where 

A=l\  -!>„]•  *=[%]•  C  =  IC»'CS1-  (8) 

For  the  elastic  dynamics  of  space  structures,  there  is  always  some  (possibly  small)  damping 
D0  appearing  in  (1)  which  causes  A  to  be  dissipative.  Thus  the  criteria  of  the  Hille-Yoshida- 
Phillips  theorem  [7]  are  satisfied  and  A  generates  a  Co-semigroup  with  an  operator  which 
we  write  suggestively  as  cAt.  Moreover,  such  models  are  ‘hyperbolic’  [9]  in  the  sense  that 
the  semigroup  is  a  contraction,  i.e.,  || e"4*  ||  <  1  and  all  but  the  zero  frequency  poles  are 
only  slightly  damped;  ||eA,j|  <  e~il  for  Borne  small  6  >  0.  We  remark  that  some  popular 
models  for  structural  elements  such  as  beams  with  material  damping  may  not  fit  in  this 
framework  [10,82].  However,  this  framework  includes  models  appropriate  for  considerations 
of  wave-like  dynamics  which  propagate  causally  in  the  spatial  domain.  For  such  models, 
the  question  of  how  to  compute  controls  u(t)  and  system  response  x(t)  focuses  on  the 
so-called  mild  solution  of  (7); 
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For  transient  disturbances,  the  standard  vibration  control  problem  is  essentially  an  op¬ 
timal  regulator  problem  which  is  readily  solved  by  a  linear  state  feedback  which  minimizes 
the  performance  index  (parametrized  by  a  real  scalor  £  >  0) 

J(u)  =  ffllvll’  +  •MV  (10) 

where  ||  •  ||  is  the  Euclidean  norm  on  the  appropriate  finite  dimensional  space.  This  is  the 
generic  control  problem  surveyed  in  Balas  [11].  In  this  paper,  we  will  concentrate  on  the 
construction  of  state  space  models  and  computational  aspects  of  equations  of  the  form  (1) 
and  of  optimal  (discrete)  controls  u,  (t)  appearing  in  (3). 

Various  methods  are  available  for  approximation  of  the  system  (7)  [12,13].  One  popular 
method  is  based  on  a  modal  (eigen-)  expansion  of  A  which  generates  a  sequence  of  finite 
dimensional  subspaces  Mk  C  U,k  =  1,2,...,  where  Ut  —  span^-,  j  =  l,...,fc}  and  the 
<i>j(z)  are  eigenfunctions  (or  mode  shapes)  for  A  [10,14].  Based  on  this  approximation,  a 
sequence  of  finite  dimensional  models  for  (5)  can  be  generated; 

*<*>(«)  =  A(fc,xW(t)  +  (11) 

Using  a  truncated  model  (11)  with  k  finite  and  the  performance  index  J(u)  projected 
onto  the  (finite  dimensional)  space  Mt,  one  can  solve  the  associated  optimal  control  problem 
for  the  first  k  modes  of  A.  This  is  a  classical  approach  and  encompasses  Ritz-Galerkin 
methods  [10]  as  well  as  spline  methods  [12,13].  However,  as  noted  in  Balas  [ll]  in  all 
but  a  few  special  cases,  the  control  law  when  applied  to  the  system  (5)  will  excite  higher 
order  modes.  The  inherent  robustness  and  stability  properties  as  well  as  the  degree  of 
suboptimality  of  control  laws  based  on  such  truncated  modal  approximations  has  received 
a  great  deal  of  attention  in  both  the  engineering  and  mathematics  literature  [9,11].  Various 
alternate  approaches  are  available  which  deal  directly  with  infinite  dimensional  control 
problem  given  by  (10)  and  (9) — at  least  abstractly  (cf.  Russell  [15]).  One  method  suggested 
by  Davis  [16,17]  offers  the  advantage  of  a  computational  procedure  for  approximating  the 
true  optimal  control  in  terms  of  the  required  control  bandwidth.  The  method  is  based  on 
an  extension  of  a  Weiner-Hopf  solution  [17]  for  the  abstract  control  problem. 

1.1.2  Wiener-Hopf  Control 

In  the  context  of  the  regulator  control  problem  given  by  the  minimization  of  (10)  subject 
to  the  infinite  dimensional  model  (7),  spectral  factorization  can  be  seen  to  provide  an 
alternate  solution  [18]  to  a  Riccati  operator  equation.  In  a  series  of  papers,  J.  Davis  and 
his  students  [16,17]  have  explored  the  application  of  spectral  factorization  for  control  design 
of  a  class  of  distributed  parameter  models  for  long  trains  with  multiple  locomotives.  Also 
in  Avramovic,  et  al.  [19]  considerations  for  application  of  these  results  to  flexible  structures 
were  given.  The  details  of  this  method  are  discussed  in  Section  2;  however,  in  this  section, 
we  will  summarize  the  relevant  results  and  highlight  the  significance  for  modeling  of  flexible 
structures. 

Taking  Laplace  transforms  in  (7)  allows  one  to  write  (at  least  formally) 

Y  (s,  z)  =  CR(a;  A)i(0,  z)  +  CR(a;  A)BU(«).  (12) 


•  F* 


br ^ a  ^ ^ ^ ^ ,  •*.  a **»  **■  **,  ■'»  •*»  a  * *.  •*.»*.»*.*«••* 
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The  transfer  function  is  G(«)  =  CR{s\A)B  where  R(a\A)  is  the  resolvent  operator  for  A , 
R(s;A)  :  M  — ►  p(A)  C  )/,  where  p(A)  is  the  resolvent  set  (or  compliment  to  the  spectrum 
of  A). 

The  optimal  control  law  which  minimizes  (10)  subject  to  (7)  consists  of  (linear)  state 
feedback 

u(t)  =  -B'Kx{t,z)  (13) 

where  B *  is  the  formal  adjoint  of  B  in  U.  We  remark  that  B'K  is  an  (linear)  integral 
operator  on  1/  which  can  be  computed  exactly  without  recourse  to  an  infinite  dimensional 
Riccati  operator  equation  via  the  formula  [17] 

B'K  =  r  [F*(«w)]"1G*(«w;  AjCJItiwjAjdh;,  (14) 

2?r  J- oo 

where  F(s)  is  the  unique,  causal  spectral  factor  of 

I  +  GT(-s)G(s)  =  F(s)Fr(-s).  (15) 

For  the  infinite  dimensional  system  (7),  the  transfer  functions  G(s),  F(s)  are  irrational  and 
F(s)  (resp.  FT(— s))  is  analytic  for  5Re  «  >  0  (resp.  Re  a  <  0).  Computational  algorithms 
for  (15)  are  given  in  Davis  [17]  and  [42]  where  a  numerical  algorithm  is  given. 

In  this  framework,  the  question  of  modeling  of  flexible  structures  for  the  design  of 
feedback  control  for  suppression  of  (linear)  vibrations  centers  on  computation  of:  (i)  the 
irrational  transfer  function  G(s),  and  (ii)  the  resolvent  R{a\  A).  Then  spectral  factorization 
in  (15)  is  performed  using  the  Davis  algorithm  [17]  and  a  numerical  approximation  to  (14) 
can  be  obtained  which  is  valid  in  the  frequency  bandwidth  of  the  desired  control  action. 
We  remark  that  although  all  computations  are  in  the  frequency  domain,  the  resulting 
control  law  is  a  linear,  distributed  state  feedback. 

In  this  report  we  start  in  Section  2  with  a  detailed  discussion  of  Wiener-Hopf  method 
for  computing  distributed  state  feedback  control.  A  numerical  algorithm  is  presented 
and  practical  aspects  of  the  method  are  considered.  In  Section  3  we  provide  a  detailed 
discussion  of  frequency  response  calculations  and  modeling  of  hybrid  structures  consisting 
of  distributed  and  lumped  elements.  In  many  cases  of  practical  interest  large  repetitive 
truss  structures  can  be  effectively  modeled  as  continua.  Section  4  provides  details  of  a 
method  called  homogenization  and  its  application  to  this  issue.  New  results  are  provided 
on  the  problem  of  control  design  based  on  the  homogenized  effective  continuum  model. 
Finally  in  Section  5  several  examples  of  modeling  and  control  computations  are  given. 
Appendices  include  computer  code  generated  as  part  of  this  project  and  a  short  discussion 
of  algorithms. 


■  ■ '  #  '  »  ’•  v  ’«*  V 
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2  Wiener-Hopf  Methods  for  Computation  of  Opti¬ 
mal,  State-Feedback  Control 

The  connections  between  least  squares  optimization,  spectral  factorization,  and  algebraic 
Riccati  equations  have  been  considered  important  in  control  theory  for  many  years.  (See 
e.g.,  Anderson  [20],  Brockett  [18],  Willems  [21],  Helton  [22],  and  the  references  therein). 
To  see  how  the  connection  arises,  consider  the  standard,  finite-dimensional,  infinite  time, 
quadratic  regulator  (LQR)  problem: 


ii«wii’+ii»wir<« 


subject  to  the  linear,  time-invariant  system  model 

x(t )  =  Ax[t)  +  Bu[t),  x(0)  =  Xo  (17) 

and  controlled  output 

y(t)  =  Cx(t),  t  >  0.  (18) 

The  transfer  function  relating  the  Laplace  transform  of  the  input  vector  u(s)  to  the  output 
vector  y(s)  is  G(s)  =  C[sl  —  A]~XB.  The  optimal  control  is  known  to  be  a  linear  state 
feedback  u(f)  =  —  Kgptxlt)  =  —BTPx(t)  where  P  is  the  unique,  positive  definite  symmetric 
solution  to  algebraic  (matrix)  Riccati  equation, 

PA  +  AtP  -  PBtBP  +  CTC  =  0.  (19) 

Standard  algebraic  manipulations  based  on  (16)-(19)  provide  the  spectral  factorization 
relation 

[/  +  Kopt{~sI-A)-lB\T[I  +  Kopt{sI-A)-lB\ 

=  /  +  Bt(-s  -  AT)~lCTC[sI  -  A)~XB.  (20) 

[23,  pp.  68]  which  we  rewrite 

H(s)  =  I  +  Gr(-s)G(s)  =  Fr(-s)F(s)  (21) 

where  F(s)  =  I  +  Kopt\sI  —  A\~XB  is  the  causal  spectral  factor  of  if(s). 

An  explicit  integral  equation  can  be  derived  for  the  optimal,  state  feedback  under  the 
additional  assumption  that  the  spectrum  of  the  operator  A  in  (17)  is  contained  the  open, 
left  half-plane  (C_). 

Theorem  1  Given  the  optimal  linear  regulator  problem  as  defined  in  (16)-(18)  the  opti¬ 
mal  feedback  gain  (if  it  exists)  can  be  computed  as 

K„t  =  ^~  l”  [F*(iw)]-1  G*(iw)CR(iutA)du>.  (22) 

2n  J — oo 

under  the  assumption  that  A  is  a  stable  (matrix)  operator. 


SEI-TR-86-14 


8 


In  (22)  we  use  the  notation  F*(iu)  =  Fr(— iu)  and  *(uj;  A)  =  [*u;/  —  A]-1  the  (matrix) 
resolvent  for  A.  One  paradigm  for  computing  LQR  type  control  laws  for  distributed 
parameter  systems  is  based  on  extending  this  computation  (rather  than  the  usual  approach 
of  finite  element  modeling  followed  by  matrix  computations).  This  requires  the  spectral 
factorization  of  H(s)  which  may  be  rational  even  though  the  system  is  infinite  dimensional. 
This  will  happen  typically  when  “modal”  control  is  the  objective  defined  in  (16);  i.e.,  the 
operator  C  maps  to  a  finite  dimensional  modal  subspace  (cf.  Section  2.3  for  details).  The 
resolvent  arising  in  (22)  can  be  computed  from  a  Green’s  function  for  the  problem.  In 
Section  3  we  provide  more  details. 

Proof  of  Theorem:  From  standard  results  [18]  the  optimal  control  (if  it  exists)  will 
stabilize  the  closed  loop  system  so  that  a  {A  —  BKopt)  C  C_  where  G_  is  the  open  left  half 
of  complex  plane.  Construct  a  closed  rectifiable  contour  T  in  the  complex  plane  consisting 
of  a  relatively  large  portion  of  the  tcj-axis  and  a  semicircular  portion  in  the  left  half  plane 
such  that  T  encircles  a  (A  -  BK^t)  in  the  positive  sense.  Let  Ae  =  A  -  BK^t.  Then 

(23) 

By  assumption  <r(— Ar)  is  contained  in  C+  and 

From  (19)  we  write 

AtP  +  PAt  =  -CTC. 

This  leads  to  the  relation 


P(sl  -  Ac)'1  +  (-sJ  -  AT)~lP  =  (-sJ  -  AT)-'CTC{aI  -  Ac)’1.  (26) 

Now  integrate  (26)  on  T  and  use  (23)  and  (24)  to  get 

P  =  -  *r)-lCTC(sI  -  A,)~'ds.  (27) 

Since  the  optimal  gain  is  K^t  =  BTP  =  \PB}T  we  get 

Kop,  =  2 hi  BT^Sl  ~  Ac)"1C,TC(-sJ  -  AT)~xd8.  (28) 

Now  from  (20)  and  (21)  we  get  that 

C(8l  -  A*)-1*  =  C(8l  -  A)"1*  [/  +  Kor^I  -  A)"1*]"1  (29) 

and  therefore 

C(s/-Ae)-1*  =  G(s)F"1(s). 

Thus  we  can  write 

=  ^/rf‘rWGTWC(-<l  -  A)'1*.  (30) 


Finally  (22)  is  determined  by  substituting  s  —  tw  under  the  observation  that  G(iu>)  — *  0 
as  w  — *  oo. 


(24) 

(25) 


e 
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For  (22)  to  hold  when  G(s)  is  irrational  we  must  consider  some  limitation  of  the  spectral 
properties  of  the  associated  infinite-dimensional  operator  A.  Consideration  for  the  compu¬ 
tation  of  the  Riccati  operator  P  by  the  integral  formula  (22)  suggests  that  the  spectrum  of 
A  must  consist  of  a  countably  infinite  set  of  eigenvalues  so  that  certain  path  integrals  can 
be  computed.  The  spectra  o(Ae)  and  o(A)  must  be  separated  by  the  imaginary  axis  (in¬ 
cluding  the  point  at  infinity).  Also,  the  observation  that  G(iu>)  -*  0  as  w  — *  oo  for  rational 
functions  must  be  replaced  by  a  similar  assumption  which  further  restricts  the  class  of  ir¬ 
rational  transfer  functions  [24].  As  it  turns  out  the  usual  assumptions  used  in  constructing 
continuum  models  for  mechanical  structures  serve  to  restrict  the  resulting  transfer  func¬ 
tions  appropriately  [51,27].  However,  such  transfer  functions  cannot  be  computed  reliably 
from  finite  element  models. 

Effectively  these  additional  assumptions  will  be  guaranteed  by  the  class  of  transfer 
functions  for  which  spectral  factorization  can  be  performed.  We  will  consider  spectral  fac¬ 
torization  for  distributed  systems  next.  More  importantly  with  the  usual  assumptions  used 
in  constructing  continuum  models  for  mechanical  structures  it  appears  that  the  resulting 
transfer  functions  will  have  the  appropriate  properties.  However,  in  this  study  we  have 
encountered  much  confusion  in  the  literature.  As  discussed  in  the  introduction  a  major 
goal  of  this  study  was  to  provide  a  consistent  method  for  model  construction  leading  to 
an  appropriate  class  of  transfer  functions  for  models  which  are  both  well-posed  and  have 
appropriate  spectral  properties.  In  the  next  section  we  delineate  the  specific  class  of  trans¬ 
fer  functions  for  which  spectral  factorization  can  be  computed  efficiently  by  a  numerical 
algorithm.  We  indicate  the  basis  for  the  algorithm  and  discuss  the  implications  of  sam¬ 
pling  and  interpolation  of  the  spectral  factor.  Finally  we  discuss  the  relationship  between 
rational  approximation  and  modal  control. 


« 


2.1  Spectral  Factorization  by  Frequency  Sampling 
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In  this  section  we  review  the  basis  for  an  interactive  algorithm  for  computation  of  a  fre¬ 
quency  sampled  spectral  factor.  The  algorithm  (due  to  Davis  and  Dickinson  [17])  provides 
an  effective  computational  tool  for  obtaining  the  optimal  gain  Kopt  via  (22)  without  regard 
to  computational  difficulties  associated  with  large  dimensional  Riccati  equations.  Conver¬ 
gence  of  the  algorithm  depends  on  certain  technical  assumptions  which  delineate  the  class 
of  transfer  functions.  In  the  finite  dimensional  setting,  a  recursive  algorithm  for  compu¬ 
tation  of  the  causal  spectral  factor  F(s)  follows  from  a  Newton-Raphson  iteration  for  the 
matrix  Riccati  equation  (19)  given  an  initial  stabilizing  feedback  Jf0  =  B*P0\ 

Pn+l{A  -  BB*Pn)  +  (A  -  BB'PnyPn+l  =  -C*C  -  PnBB*Pn. 

At  the  nth  iteration  one  can  take  an  approximation  to  the  causal  spectral  factor  as 

Fn(s)  =  /  +  H*P„(s/-A)-,H. 

Then  following  Davis  and  Dickinson  [17]  this  leads  to  the  form  of  the  algorithm 

=  pt  (3i) 
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where  P+  is  the  causal  projection  operator  defined  on  the  convolution  algebra  J  ©  £1  or 
/  ©  £  2  by 

P+  {  J  +  /“  /(t)e-<w‘dt}  =  /  +  jH 

Computation  of  causal  projection  of  a  signal  is  a  standard  problem  in  signal  processing 
which  can  be  effectively  solved  through  the  use  of  a  Hilbert  transform  [25],  [26]. 

Definition  2  The  Hilbert  transform  of  a  signal  f(t)  is  given  as  [25] 


The  basic  utility  of  the  Hilbert  transform  is  evident  by  examination  of  its  Fourier  transform. 
Let  F(u>)  be  the  Fourier  transform  of  7(f)  and  F(u)  be  the  Fourier  transform  of  /(f). 


P(u>)  =  [°°  -  [°°  l^-dr  t^'di. 
J- 00  *  J — oo  t  —  T 


Since  for  Z{u)  —  ~e,u,tdt  we  see  that  H(u)  is  a  complex  function  with  the  properties 


|H(u>)|  =  1 

arg H(u)  =  I  -'/*•  w>°  . 

*  1  '  (  x/2,  u>  <  0 


Therefore  from  (33)  and  (34)  we  get  that 


p(u)  =  J  w  >  0 

'  *  \  iF(w),  u  <  0 


Now  to  compute  the  causal  projection  of  /(f)  (as  in  (31))  given  frequency  domain  data 
F(u)  we  first  compute  the  Hilbert  transform 

Flu)  =  F(w)  •  — . 

rui 

By  duality  of  the  Fourier  transform  pair  the  property  (35)  holds  in  the  time  domain  f; 

7~x  {^(w)}  =  -i  sgn(f)/(t). 

Finally  causal  projection  can  be  computed  as 

K  {FM}  =  \  [fH  +  .*(«)]  .  (36) 

Before  we  consider  computational  issues  further  we  review  the  extension  of  this  algorithm 
to  irrational  transfer  functions.  The  questions  of  existence  and  uniqueness  of  the  spectral 
factorization  of  the  transform  H{s)  =  /  +  Gr(s)G(s)  naturally  lead  to  conditions  for  which 
5(w)  =  GrT(— »u>)G(tw)  is  positive  semidefinite  for  u  real.  In  the  convergence  proof  of  the 
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iteration  (31)  Davis  assumes  that  G(s)  is  the  Fourier  transform  of  a  real,  vector-valued 
function  which  is  both  integrable  and  square  integrable;  i.e.,  G(iu>)  €  J{£1  D  £2). 

With  these  assumptions  it  is  clear  from  the  classical  theory  of  Gohberg  and  Krein  [28] 
that  tf(s)  has  a  unique  spectral  factorization  as  given  in  (21)  with 

T{£+) 

where  J(£j")  is  the  class  of  Fourier  transforms  of  functions  in  £\  with  positive  support. 
As  noted  in  [17],  the  assumptions  on  G(a)  in  fact  imply  that 

F±1(iu;)-/e7(£J-n£+) 

and  F(iuj)  =  F(—iu). 

Using  this  assumption  Davis  is  able  to  show  that  the  recursion  (31)  starting  from  an 
initial  F0(M  €  T{£2  Cl  £2)  has  all  iterates  F„(*u>)  G  J(£j  ft  £2)  and  that: 

lim^oo  F*(u>)Fn{u)  =  H(u)  almost  everywhere  (37) 

limn^oo  F„(s)  =  F(s)  for  all  s  with  Die  s  >  0.  (38) 

These  results  should  not  be  surprising  for  the  class  of  transfer  functions  considered.  The 
following  theorems  summarize  well  known  properties  of  £\  and  £2  functions. 

Theorem  3  ([29])  If  /  e  £\  then 

1.  w>-*  /(u>)  is  uniformly  continuous  for  w  6  Di 

2.  /M  e  £«, 

3.  I/Ml  — »  0  as  |u;|  — ♦  oo 

4.  f(t)  =  A  f(iu})e'utdu}  almost  everywhere  in  Di. 

Theorem  4  (Parseval’s  theorem  [29])  If  /  £  £j  then 
1. 

/_”  I/ll*  =  £  /“  l/MI<iw 

2.  as  N  — ►  oo 

„  m'-"*  -ts.  niu) 

3.  as  N  — »  oo 

h  £  m 

So  we  see  that  /  6  £\  means  that  /(»w)  is  bounded  on  u /  and  therefore  the  Fourier 
transform  is  well  defined  while  /  G  £2  provides  consistent  approximation  theory  for  band 
limited  signals.  We  note  that  the  third  £\  property  means  that  such  transfer  functions  are 
effectively  band  limited  and  strictly  proper. 


O'  a."  a"  a*  a*  a*  a*  a  *  ».  *  a  *  .  •  .  i  . 
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Application  of  spectral  factor  to  the  integral  formula  (22)  will  further  restrict  consid¬ 
eration  to  G(a)  both  causal  and  stable  so  that  G(«)  is  analytic  for  Re  a  >  0.  Thus  the 
transfer  functions  we  are  considering  belong  to  the  Hardy  space  G(s)  €  H2  n  H°°.  Recall 
that  by  definition  /  £  H2  if  /  is  complex  valued  and  analytic  in  C+  (the  open  right  half 
of  the  complex  plane)  and 

roe  A 

sup  /  | /(<r  +  iu)\du>  <  oo. 

o>0  J— oo 

while  /  £  H°°  if  /  is  complex  valued  and  analytic  in  C+  and 

sup  |7(s)|  <  00. 

•ec+ 

The  first  property  is  inherited  from  f  €.  £  j  while  the  second  comes  from  f  €  £ j. 

Finally,  we  remark  that  for  mechanical  structures  that  transfer  function  models  can  be 
factorized  in  a  “product  expansion”  [30] 

>ia\  =  TT  1  +  (*V*t)  (39\ 

/l)  Mi +(•*/«»■  (39) 

Thus  we  conclude  these  remarks  by  indicating  that  the  class  of  transfer  function  models 
for  which  spectral  factorization  can  be  computed  by  the  present  method  are  meromorphic 
and  have  representations  as  (39). 

2.1.1  Remarks  on  Algorithm  Construction 

With  regard  to  the  integral  formula  for  the  optimal  state  feedback  gain  it  is  clear  that  we 
need  [jF(iw)]-1.  Thus  as  suggested  by  Davis  it  is  convenient  to  implement  the  iteration  in 
the  form 

=  [F„r  (/ + k  (ra-  »  mr1  -  /})"’ .  (<o> 

Furthermore  by  initializing  with  Fo  a  diagonal  matrix  with  diagonal  elements  equal  to  the 
spectral  factors  of  the  diagonal  elements  of  H  the  second  term  of  (40)  remains  a  pertur¬ 
bation  of  the  identity  (since  [F*]-1ff[F„]_1  —  /  — ♦  0)  which  regularizes  the  computations. 
The  diagonal  initialization  guarantees  that  the  first  residual  [FoJ'^lFo]-1  has  ones  on  the 
diagonal  and  all  off  diagonal  elements  less  than  one  in  magnitude.  Using  the  properties  of 
the  Hilbert  transform  and  the  formula  (36)  one  can  readily  compute  the  causal  spectral 
factor  for  the  individual  scalor  transfer  functions  directly  (i.e.  without  iteration).  In  par¬ 
ticular,  the  ktK  diagonal  element  of  /f(u;),h*(u>)  is  a  real  valued  function  with  hk(u)  >  0. 
Let  hk(u)  be  the  Hilbert  transform  of  hfc(w);  viz., 


M»)  -  i  r 

K  J  —oo  U)  —  O 

then  the  causal  spectral  factor  h*(iu>)  =  /*(— »"w)/*(»w)  is  given  by 

A(w)  =  v/ReMw)«"lhM. 


r.  •  .  *  .  •  «  4  •  *  *  .  *  •'  •  •  •  .*  #  v"  v*  V  s'  • 
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Finally,  we  remark  that  numerical  computation  of  the  Hiblert  transform  can  be  achieved 
in  several  ways.  Direct  numerical  integration  of  (41)  is  complicated  by  the  fact  that  the 
integral  is  convergent  in  the  Cauchy  principal  value  sense.  Effective  quadrature  algorithms 
for  such  problems  have  been  coded  and  tested.  A  public  domain  version  utilizing  an 
adaptive  quadrature  procedure  is  available  in  the  routine  QAWC  contained  in  a  software 
package  called  QUADPACK  (31]. 

Another  approach  is  taken  by  Davis  [17]  who  uses  an  algorithm  of  F.  Stenger  [32].  This 
procedure  essentially  implements  a  discrete  (sampled)  version  of  the  required  computation 
using  a  digital  Hilbert  transform.  For  a  finite  number  of  sample  points  the  Hilbert  trans¬ 
form  computation  can  be  implemented  by  taking  a  discrete  “fast  fourier  transform”  (FFT) 
of  the  sampled  data  and  shifting  the  imaginary  part  of  the  transformed  data  according 
to  (33).  It  is  well  known  that  control  of  error  induced  by  the  sampling  process  (Gibbs 
phenomenon)  requires  the  careful  choice  of  “data  windows”  or  weighting  functions  for  the 
computation.  Although  these  considerations  are  well  described  in  the  literature  on  digital 
signal  processing  [26],  [34],  [33]  it  is  not  apparent  that  these  considerations  were  contained 
in  the  work  of  F.  Stenger. 

In  our  experience  the  direct  implementation  of  the  algorithm  described  in  Davis  and 
Dickinson,  based  on  the  causal  projection  of  F.  Stenger  [32]  is  unreliable  at  best.  Since 
detailed  design  of  window  functions  for  discrete  Hilbert  transforms  was  considered  outside 
the  scope  of  the  present  study  we  have  found  it  expedient  to  use  the  adaptive  quadrature 
software  from  QUADPACK  [31].  It  should  be  noted  however  that  this  approach  may 
suffer  in  production  grade  applications  by  the  computational  inefficiency  of  the  quadrature 
procedure  by  comparison  with  an  FFT  implementation  of  the  discrete  Hilbert  transform. 


2.2  Importance  of  Damping  in  Models  for  Distributed  Control 
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For  the  present  study  we  have  attempted  to  compute  distributed  controls  for  several  ex¬ 
amples.  Several  negative  results  lead  us  to  reevaluate  the  theoretical  basis  for  control 
and  modeling  problem.  Our  difficulties  stem  from  precise  computation  of  the  irrational 
transfer  functions  for  several  distributed  models  which  will  be  discussed  later.  Thus  our 
control  laws  are  computed  and  tested  on  bona  fide  distributed  parameter  models  and  not 
on  finite  dimensional  approximations;  therefore  the  idiosyncracies  of  many  available  mod¬ 
els  (which  we  unfortunately  tried  first)  lead  to  difficulties.  In  this  section  we  briefly  review 
our  most  recent  conclusions  about  the  assumptions  in  modeling  for  mechanical  structures 
and  resulting  computation  of  optimal  control  laws. 

The  most  pertinent  comments  in  the  literature  appear  to  be  summarized  in  the  work 
of  Gibson  [10].  In  particular,  the  importance  of  damping  and  the  ability  to  compute 
distributed  control  laws  for  infinite  dimensional  models  for  mechanical  structures.  It  is 
true  that  computation  of  control  laws  for  such  problems  will  always  involve  approximation 
in  modeling,  control  problem  formulation,  and  numerical  computation.  Gibson  focuses  on 
this  question  of  approximation  in  modeling  and  control  for  distributed  parameter  problems 
via  the  use  of  finite  dimensional  models  and  control  computation  based  on  these  models. 
Although,  on  the  surface,  the  Davis  method  (based  on  spectral  factorization  of  irrational 
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transfer  functions)  is  not  subject  to  finite  dimensional  approximation  (at  least  in  the  sense 
of  modal  expansions)  it  is  clear  that  approximation  via  sampling  and  interpolation  of  the 
transfer  functions  involved  is  required.  The  examples  in  Section  5  suggest  that  additionally 
some  form  of  rational  approximation  may  also  be  required  for  numerical  computation  of 
the  distributed  gain. 

Let  us  consider  some  observations  of  Gibson  [10].  First,  Consider  the  class  of  elastic 
mechanical  structures  modeled, 

x(t,z)  +  C0i(t,z)  -f  A0x(t,z)  =  Bu(t)  (43) 

where  x  6  V.  (H)  appropriately  chosen  to  match  required  boundary  conditions,  u(t)  is  a 
finite  dimensional  vector  of  controls,  A<>  is  a  self  adjoint  operator  Ao  :  D(Ao)  — ►  "H  densely 
defined  on  M.  Gibson  further  assumes  that  Ao  is  coercive ;  i.e. 

(Ao*,*)#  >  p’NIi,  x  £  D(A0) 

and  Aq1  is  compact.  C0  is  nonnegative,  symmetric  linear  operator  and  there  exists  -y  >  0 
such  that 

\\C0x\\  ^  ^llAox]]#,  X  e  D(Ao).  (44) 

Finally,  Bo  is  taken  to  be  a  bounded  operator.  Gibson  shows  that  (44)  is  a  necessary 
condition  for  the  resulting  semigroup  operator  for 


[  -Ao  -Co  J 

on  H  x  )i  to  be  uniformly  exponentially  stable;  i.e.,  there  exists  M  >  0,a  >  0  such 
that  ||eA‘||  <  Me~at  for  t  >  0.  Such  behavior  corresponds  to  the  case  when  damping 
provides  a  uniform  decay  rate.  Gibson  addresses  the  question  of  convergence  and  stability 
of  a  sequence  of  finite  dimensional  (possibly  modal)  approximate  control  problems  to  the 
unique  optimal  control  for  the  distributed  model.  His  results  indicate  [10,  theorem  4.1] 
that  for  the  quadratic  linear  regulator  problem  for  the  distributed  model  above  (43)  that 
if  a  solution  exists  the  exact  optimal  feedback  control  provides  a  closed  loop  system  whose 
infinitesimal  generator  Ad  =  A  +  BK^t  generates  a  strongly  continuous  semigroup  which 
is  uniformly  exponentially  stable. 

For  systems  without  damping;  i.e.  C0  =  0  in  (43),  Gibson  shows  that  there  can  be 
no  nonnegative,  selfadjoint  solution  for  the  algebraic  Riccati  equation  for  the  distributed 
model.  This  follows  from  the  observation  [10,  theorem  5.13]  that  for  Co  =  0,  the  semigroup 
generated  by  A  +  £  for  any  £  a  compact  linear  operator,  cannot  be  uniformly  exponentially 
stable. 

For  systems  with  damping  Gibson  shows  that  the  damping  must  be  such  that  A  (the 
generator  for  the  open  loop  system)  is  uniformly  exponentially  stable  in  order  that  uniform 
exponential  stability  of  the  closed-loop  system  can  be  obtained  by  compact  linear  feedback. 
The  significance  of  this  fact  is  that  the  physical  nature  of  damping  must  be  considered  in 
the  computation  of  distributed  parameter  control  and  for  problems  where  the  available 
control  affectors  provide  ‘localized’  (in  the  spatial  domain)  forces  (or  torques).  Damping 
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with  (at  least)  a  uniform  decay  rate  (for  essentially  all  modes)  is  required  for  such  systems 
to  be  controlled  with  uniform  exponential  stability.  Although  it  is  generally  agreed  that 
physical  structures  will  have  this  property  it  is  somewhat  disconcerting  that  many  standard 
models  for  simple  structural  elements  like  beams  with  internal  damping  do  not  have  these 
properties.  In  fact,  with  the  exception  of  viscous  damping  there  does  not  appear  to  be  a 
single,  universally  accepted,  well-defined  mathematical  model  for  beam  dynamics  which  is 
obviously  appropriate  for  the  study  of  distributed  parameter  control  of  structures  [82]. 

In  Section  3  we  will  discuss  damping  models  for  structural  elements  in  detail.  At  this 
point  we  remark  that  in  this  study  it  became  apparent  that  the  performance  of  active 
control  of  elastic  structures  was  seen  to  be  heavily  dependent  on.  the  type  of  damping 
models  used.  If  one  is  willing  to  assume  a  sufficient  amount  of  viscous  damping  then 
stabilizing  control  can  be  computed.  Since  viscous  effects  will  definitely  not  be  present  in 
space  applications  one  is  faced  with  a  choice  of  several  models  for  internal  damping.  Many 
of  these  models  lead  to  problems  for  which  either  a  uniform  decay  rate  is  not  available 
or  for  which  the  spectrum  of  the  operator  contains  more  than  mere  eigenvalues.  These 
models  cannot  be  stabilized  by  compact  linear  feedback  whether  it  is  computed  by  Weiner- 
Hopf  methods  or  by  modal  approximation  and  solution  of  a  finite-dimensional  quadratic 
regulator  for  the  reduced  model.  However  if  one  takes  the  second  path  the  relevant  stability 
questions  are  completely  lost  in  the  model  reduction. 

What  is  an  appropriate  choice  for  modeling  internal  damping  in  space  structures  i8 
definitely  an  open  question.  This  issue  has  apparently  not  received  a  great  deal  of  attention 
in  the  aerospace  industry  especially  in  that  portion  of  the  community  involved  in  control 
of  large  flexible  structures.  The  standard  procedure  here  (as  evidenced  in  for  instance  the 
ACOSS  program)  is  to  assume  “low”  modal  damping  can  be  added  to  the  reduced-order 
model  prior  to  control  system  design.  In  the  present  frequency  domain  computations  one 
is  forced  to  resolve  these  issues  before  proceeding  with  control  computation. 

For  the  computation  of  distributed  control  according  to  (30)  it  is  therefore  necessary 
that  damping  be  of  (at  least)  uniform  exponential  type.  However,  even  more  is  required 
in  practice.  Since  Kopt  will  be  computed  (in  our  approach)  by  numerical  evaluation  of 
(30)  it  is  clear  that  the  integration  will  be  performed  over  a  finite  bandwidth  by  numer¬ 
ical  quadrature.  The  required  bandwidth  for  integration  depends  crucially  on  the  open 
loop  system  bandwidth  for  the  continuum  structural  model.  This  means  that  effectively 
damping  must  be  ultimately  proportional  to  frequency  for  all  high  frequency  modes.  This 
is  consistent  with  the  usual  notion  of  material  damping  but  is  required  here  for  numerical 


reasons. 


2.3  Practical  Aspects  of  Computing  Distributed  Control  Laws 

In  this  section  we  review  the  computational  aspects  of  Wiener-Hopf  control  based  on  spec¬ 
tral  factorization  and  distributed  gain  computation  and  its  relation  to  standard  methods 
based  on  finite  dimensional  approximation  of  the  evolution  equations.  We  highlight  var¬ 
ious  methods  for  model  approximation  (i.e.  reduced  order  modeling).  The  methods  we 
employ  are  based  on  frequency  domain  computations  for  state  space  models  and  can  be 
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compared  directly  with  state  space  methods.  However,  as  discussed  in  Youla  et  al  [35]  the 
Wiener-hopf  method  can  treat  models  which  have  no  state  space  analog  [24]. 

Reduced  order  modeling  (ROM)  of  the  abstract  control  problem  (7)-(10)  can  proceed 
in  many  ways  [36].  At  a  fundamental  level  analysis  of  ROM  involves  decomposition  of 
the  state  space  ^(O)  via  an  orthogonal  decomposition  M  =  Ms  ©  Mr  where  Mn  is  a  finite 
dimensional  subspace  of  1/(0)  and  Mr  is  a  residual  space.  For  example,  with  x(z)  (E  #(0) 
let  x  be  expanded  as 

*(*)  =  f 

t=i 

which  we  assume  converges.  The  series  {<f>k{z))T=i  provides  an  orthogonal  basis  for  1/(0). 
Then  the  ROM  space  is 

K,v(n)  =  span{<£i(z),$2(2),. . .  ,4>n{z)}- 

This  approach  embodies  finite  element  methods  based  on  Ritz-Galerkin  approximation  as 
well  as  truncated  modal  approximation. 

In  Balas  [36]  tradeoffs  are  discussed  for  the  formulation  of  optimal  control  for  dis¬ 
tributed  parameter  systems  (DPS)  by  projection  onto  reduced-order  6ubspaces.  In  Bern¬ 
stein  and  Hyland  [37]  an  optimal  projection  approach  is  described  for  reduced-order  models. 
In  applications  it  may  be  desirable  to  define  a  control  objective  of  the  form  (10)  so  that 
only  a  finite  number  of  modes  are  considered  for  active  control.  Typically  the  accuracy 
of  any  mathematical  model  degrades  with  higher  frequencies.  Sensors  and  actuators  have 
finite  bandwidth. 

Let  PN  be  an  orthogonal  projection  Pn  :  1/(0)  — ♦  1/^(0).  The  quadratic  control 
objective  (10)  is  defined  with  respect  to  a  controlled  output  y(t)  =  Cx(t,z)  given  by  the 
operator  C  :  1/  — ►  5RP.  A  reduced-order  control  objective  can  be  defined  by  the  choice 
C  —  Pjv.  For  the  finite  dimensional  version  of  the  control  problem  (7)-(l0)  it  is  known  that 
an  optimal  control  law  exists  which  stabilizes  the  system  if  and  only  if  (A,  B )  is  stabilizable 
(i.e.,  any  state  trajectory  contained  in  the  uncontrollable  subspace  is  exponentially  stable) 
and  (A,C)  is  detectable  (i.e.,  any  state  trajectory  contained  in  the  unobservable  subspace 
is  exponentially  stable.)  For  distributed  parameter  models  for  elastic  structures  we  assume 
the  operator  A  has  discrete  spectrum  and  is  at  least,  uniformly  exponentially  stable  for 
essentially  all  modes.  Thus  the  optimal  control  problem  (7)— (10)  with  controlled  output 
given  by  C  =  Pn  has  a  unique  stablizing  control  law  given  by  (13). 

We  remark  that  the  formulation  is  technically  distinct  from  the  methods  discussed 
in  [36]  in  that  the  control  problem  incorporates  ROM  as  part  of  the  control  objective 
rather  by  approximation  (ROM)  of  the  (constraint)  equation  (7)  by  Ritz-Galerkin  or  modal 
truncation.  In  terms  of  the  mild  solution  (9)  we  see  that  with  y{t)  —  P n  x(f,z), 


y(t)  =  PjveA‘x(0,z)  +  f  PNeA^~o)Bu{o)da. 

Jo 


Thus  the  control  problem  can  be  equivalently  stated: 
Problem  minimize  the  objective 


rw 

Jo 


+  Hull2  dt 
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subject  to  (45). 

The  projection  Pn  is  orthogonal  so  that  for  any  x(z)  6  M  (fi)  we  can  write  P*x  4-  Prx 
where  Pr  is  the  orthogonal  complement  projection;  i.e., 

PRx{t,z)=  52  xk{t)<f>k(z). 
k=N+ 1 

Thus  (45)  can  be  written  alternately  as  an  evolution 

y{t)  =  PNAPNx(t,z)  +  PNAPRx{t,z)  + 

Roughly  put,  the  term  PsAPRx{t,z)  does  not  contribute  to  the  control  objective  so  that 
the  control  problem  is  equivalent  to  the  problem 

Studio  MOII’  +  llvW* 

subject  to 

ir(0  =  PffAy(t)  4  PjvBu(t); 
i.e.,  a  finite-dimensional  control  problem. 

The  Wiener-Hopf  approach,  taken  in  this  study,  with  the  reduced-order  control  ob¬ 
jective  further  clarifies  that  the  underlying  control  problem  is  finite-dimensional.  One 
way  to  view  ROM  in  the  frequency  domain  is  to  consider  rational  approximation  of  the 
transcendental  functions.  Consider  the  transform  of  (9); 

x(s,  2)  =  P(s;  A)x°{z)  +  R{a\A)B  u(s). 


In  Section  3  we  show  that  for  elastic  structures  this  abstract  equation  can  be  written  in 
the  form 

x(a,z)  =  J^Gr(a\ztw)x0(w)dw  +  Jfac  («,*)«(*), 

where  Gr(s;  2,  tv)  is  a  Green’s  function.  The  resolvent  operator  R(s;  A)  is  an  integral  oper¬ 
ator  with  kernel  a  Green’s  function.  Let  uk,$k  be  the  modal  frequency  and  damping 
for  A  with  respect  to  some  ROM  basis  for  )iff{ H).  Then  (at  least  on  the  surface)  modal 
truncation  can  be  viewed  as  a  truncated  partial  fraction  expansion  (see  Wie  and  Bryson 
[30]); 

Py  1(3,2)  =  Gr,N(*\z,v>)x0{w)dw  4  HBc,n(s;  z,  w)u(s) 


where 


Gr,/f(s;z,w) 


y  tl>k[z,w) 

«2  +  2 f*o»*s  +  u>2 


Hbc,n{*\  x) 


y  (t(x,w) 
iT[  4  2 f*w*s  +  wl  ’ 
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The  computational  procedure  for  the  Wiener-Hopf  problem  outlined  in  Section  2  starts 
with  the  definition  of  the  transfer  function  G(s)  from  the  control  tt(s)  to  the  controlled 
output  y(s).  For  the  case  of  ROM  control  objective  we  get 

G{s)  =  PNR{s-tA)B 

=  Hbc,n{s\z) 

a  rational  transfer  function  in  the  complex  frequency  s.  Several  algorithms  exist  for  this 
classical  problem  [38,39].  Then  the  optimal  gain  formula  requires  integration  as 

Kopt  =  f°°  \F,{iu)}-1G‘{iu)CR{iuj-,A)doj  (46) 

27 r  J — oo 

=  J  \F‘(iu)]~1G‘(iuj)C  (^J  Gr(iu>-,  z,tv)  •  dw^jdu. 

We  remark  that  the  feedback  “gain”  is  an  integral  operator  with  kernel  we  call  the  dis¬ 
tributed  gain  Kdi,  where 

Kdt.(w)  =  ±-  r  (F'(*w)]'lG*(tw)CGr(tw;*,u;)  dw  (47) 

27T  J  —  oo 

can  be  computed  as  a  path  integral  in  the  complex  plane  and  replacing  the  Green’s  function 
Gr  with  its  ROM  Gr<N  and  G  replaced  with  Gn  in  (47)  shows  that  the  integrand  is  a  rational 
function  in  s.  Thus  Kd„(w)  can  be  evaluated  by  residues. 

The  main  point  to  be  emphasized  is  that  the  numerical  procedure  discussed  in  Section  2 
takes  a  fundamentally  different  approach  by  using  frequency  sampling  in  place  of  rational 
approximation  for  the  required  computation  of: 

•  spectral  factorization 

•  distributed  gain  computation 

The  numerical  procedure  can  provide  an  effective,  computationally  efficient  alternative — 
even  for  large  finite-dimensional  problems  which  arise  in  optimal  cont-ol  of  DPS  with 
objective  which  weight  a  ROM  subspace. 

Finally,  we  conclude  this  section  with  some  remarks  about  the  practical  aspects  of 
frequency  sampling  and  numerical  computations.  The  spectral  factorization  algorithm 
has  already  been  discussed  in  detail.  The  numerical  evaluation  of  Kdi,(w)  proceeds  by 
a  trapezoidal  quadrature  over  a  finite  bandwidth  — u>0  <  w  <  w0.  The  choice  of  u;0  for 
numerical  integration  depends  on  the  bandwidth  of  the  transfer  function  G(s).  For  the  case 
of  reduced-order  control  objective  u;0  can  be  readily  determined  by  reference  to  the  rational 
transfer  function.  For  the  irrational  G(s)  it  is  necessary  that  the  high  frequency  modes  be 
sufficiently  damped  so  that  finite  bandwidth  integration  yields  a  reasonable  approximation 
to  the  true  Kdi,(w).  The  examples  in  Section  5  serve  to  illustrate  that  in  some  cases  even 
more  seems  to  be  required  for  effective  convergence  of  numerical  algorithms.  Depending 
on  the  choice  of  the  control  objective  it  seems  that  rational  approximation  may  still  be 
required  for  the  finite  bandwidth  computation  of  the  distributed  gains. 


The  frequency  domain  approach  offers  several  alternatives  for  model  order  reduction, 
approximation,  and  control  computation.  Direct  approximation  of  transfer  functions  by 
truncation  of  product  expansions  is  discussed  in  Wie  and  Bryson  [30].  This  method  has 
also  been  used  for  simulation  [40].  Here  we  expand  an  irrational  transfer  function  (assumed 
to  be  meromorphic) 

TT  Wfc(6) 

ii  M*) 

where  nfc(s),d*(s)  are  (relatively  low  order)  polynomials  in  s.  By  truncation  of  (48)  we 
effectively  match  a  finite  number  of  both  poles  and  zeros  of  the  irrational  transfer  function. 
This  method  is  clearly  distinct  from  the  modal  and  finite  element  methods.  Approximation 
of  the  zeros  of  the  transfer  function  may  be  extremely  important  for  control  system  design 
since  performance  of  designs  based  on  poorly  approximated  zero  locations  can  be  very 
sensitive  to  model  errors  [30,24].  Zero  locations  arise  in  flexible  structure  control  problems 
from  the  location  and  type  of  actuators  and  sensors. 

Finally,  in  this  report  we  have  not  addressed  the  design  of  dynamic  control  laws.  This 
question  is  discussed  in  detail  in  [24].  We  remark  that  distributed  state  feedback  control 
may  be  quite  realistic  for  the  next  generation  of  spacecraft  with  flexible  structures.  Recent 
technology  developments  in  instrumentation  for  distributed  measurement  clearly  shows  a 
trend  to  provide  low  cost  alternatives  for  distributed  measurements  of  displacements  and 
angular  rates  as  well  as  various  physical  quantities  such  as  strain  rates  [41].  In  Section  3 
we  consider  choice  of  state  space  coordinates  for  continuum  models  of  structures  including 
the  physical  significance  of  the  state  variables. 
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3  State  Space  Models  for  Distributed  Structural  El¬ 
ements 

3.1  Standard  Forms  for  Linear  PDEs 

In  this  section,  we  will  consider  the  problem  of  deriving  standard  state-space  descriptions 
for  typical  distributed  elements.  In  Section  3.2  we  use  these  descriptions  to  compute  the 
complete,  frequency  domain  response  as  required  for  the  control  problem  outlined  above. 
It  is  our  contention  that  the  systems  of  interest  to  us — specifically  beams  with  one  space 
variable  and  perhaps  several  degrees  of  freedom — can  be  represented  by  one  of  two  stan¬ 
dard  forms.  Once  identifying  the  structure  of  these  standard  models,  it  is  straightforward, 
although  far  from  trivial,  to  mechanize  the  construction  of  the  required  transfer  matrices 
using  symbolic  computation  [42].  Moreover,  in  order  to  assemble  hybrid  system  mod¬ 
els  by  the  interconnection  of  components  or  subsystems,  it  is  essential  to  have  a  clear 
understanding  of  the  causal  requirements  of  the  component  mathematical  models.  The 
following  paragraphs  develop  the  required  concepts  in  terms  of  commonly  used  structural 
elements.  Since  typical  elements  interact  at  physical  boundaries,  our  foremost  concern  is 
with  the  formulation  of  appropriate  boundary  conditions  for  well-posed,  initial-boundary 
value  problems. 

Before  proceeding,  we  establish  some  basic  terminology  associated  with  systems  of  par¬ 
tial  differential  equations.  Consider  the  system  of  first-order,  partial  differential  equations 
defined  for  t  >  0  and  0  <  z  <  L 

dw  -  dw 

E^r =  Fj^+Hw'  weRn’  (4°) 

If  E  is  nonsingular,  then  (49)  can  be  written 

dw  „dw  __  . 

m=FT*+Hw  (50> 

where  F  =  E~lF,  H  =  E~lH .  If  F  has  only  real  eigenvalues  and  a  complete  set  of 
eigenvectors,  then  the  system  is  said  to  be  hyperbolic  (see,  for  example,  Zauderer  [43]). 
If  there  are  multiple  real  eigenvalues  and  less  than  a  complete  set  of  eigenvectors,  then 
the  system  is  of  (partial)  parabolic  type.  If  all  of  the  eigenvalues  are  complex,  the  system 
is  of  elliptic  type.  Systems  with  complex  eigenvalues  are  not  causal.  Lyczkowski,  et  al. 
[44],  and  Sursock  [45]  provide  an  interesting  discussion  of  this  point  in  connection  with  a 
fluid  flow  problem.  The  underlying  problem  is  that  systems  with  complex  eigenvalues  are 
not  well-posed  as  initial  value  problems,  John  [46],  Lax  [47].  We  will  not  consider  such 
problems  any  further. 

If  E  is  singular,  (49)  can  give  rise  to  mixed  systems  of  all  types.  Some  examples  can  be 
found  in  Firedly  [48]  and  Lapidus  and  Pinder  [49].  Our  interest  in  this  case  will  be  limited 
to  purely  parabolic  systems  of  the  type 

dw  d2w  dw 

*=G7&  +  Fa;  +  H"  <»> 


.  -  * .  .  »,• 
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which  commonly  arise  in  engineering  problems. 

When  the  equations  of  motion  for  structural  elements  are  derived  from  conservation 
laws — in  particular,  from  variational  principles — the  resulting  equations  are  typically  of 
hyperbolic  type  (see,  for  example,  Crandall,  et  al.  [27]).  However,  further  standard  as¬ 
sumptions  and  approximations  reduce  the  equations  to  parabolic  systems  in  the  form  of 
(51).  In  the  following  paragraphs,  several  examples  will  be  given.  In  summary,  we  are 
primarily  interested  in  hyperbolic  systems  (49)  and  parabolic  systems  (51).  In  addition  to 
equations  (49)  or  (51),  there  are  associated  initial  and  boundary  conditions.  For  equation 
(49),  these  conditions  take  the  general  form 


initial  conditions  tu(z,0)  =  f(z) 
boundary  conditions  Eiu/(0,t)  +  Tiw(L,t)  =  g(t) 


(52) 


«v 

A 

A 


where  dim(</)=dim(u>);  and  for  equation  (51),  they  take  the  general  form 

initial  conditions  11/(2, 0)  =  f(z )  ,  . 

boundary  conditions  Eitt/(0,t)  +  £2^7(0,*)  +  Tiw[L,t)  +  rj|j(.L,t)  =  y(t)  ^  ' 

where  dim(y)=2dim(u/). 

It  is  well  known  that  the  coefficient  matrices  in  (52),  (53)  must  satisfy  certain  con¬ 
straints  if  the  problem  formulation  is  to  be  well-posed.  In  the  hyperbolic  case  (equations 
(49)  and  (52)),  these  constraints  essentially  require  that  the  boundary  conditions  be  com¬ 
patible  with  the  wave  directions.  Further  discussion  can  be  found  in  Russell  [15]  and 
Agarwala  [13].  In  the  following  sections  we  will  discuss  some  standard  models  for  beams. 
The  following  notation  is  standard  and  assumes  the  elastic  beam  is  uniform  (i.e.,  the 
parameters  are  independent  of  the  spatial  coordinate.) 


V 

v 

8 

X: 


Standard  Notation  for  Uniform  Beam  Formulae 

#cG  effective  shear  modulus 

p  mass  density 

A  cross  sectional  area 

E  modulus  of  elasticity 

I  moment  of  inertia 

L  beam  length 

z  longitudinal  coordinate 

r)(z)  lateral  displacement 

<t>  angular  rotation  of  cross  section 


5* 

V 


§ 


3.1.1  The  Timoshenko  Beam  Model 

We  will  show  how  some  conventional  beam  models  can  be  reduced  to  the  standard  formB 
described  in  the  preceding  paragraphs.  In  particular,  we  will  begin  with  the  Timoshenko 
model  and  then  consider  two  commonly  used  approximations  which  can  be  derived  from  it, 
the  Euler-Bernoulli  model  and  the  ‘string’  model.  The  equations  of  motion  can  be  derived 
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using  Lagrange’s  equations  [27]  and  in  the  absence  of  dissipation  take  the  form 


(£-*) 


along  with  the  natural  boundary  conditions  for  a  =  0,  L 

displacement  or  shear  force  .  . 

.,(«,()  =».(<)  *CA  *  > 

and 

rotation  or  moment  .  . 

El  («£!))  =  r„(t).  (56) 

The  two  equations  (54)  can  be  replaced  by  four,  first-order  equations  by  introducing 
two  new  variables,  v[z,t)  and  7(2,*): 


dt\  _  du 
dt  ~  Ite* 

i  -  »» 

^7  A 

dt  ~  dz  +  1V' 

#7  _  E  d<f> 
dt  p  dz’ 

These  equations  clearly  represent  a  hyperbolic  system  and  the  natural  boundary  conditions 
become  for  a  =  0,L 


displacement 

or  shear  force 

(58) 

q(a,t)  =  i?«(t) 

v(a,t)  =  va{t),va(t)  = 

rotation 

or  moment 

(59) 

<f>(a,t)  =  ^a(t), 

7(<M)  =  7a(0.  'MO  = 

Note  that  the  boundary  conditions  applied  to  the  first  order  system  (57)  require  the 
time  integral  of  boundary  forces  or  moments  applied  to  the  beam.  It  is  easy  to  confirm  that 
the  transfer  functions  relating  forces  or  moments  to  displacements  or  rotations  as  derived 
from  either  equations  (54)  or  (57)  are  indeed  identical  and  that  the  required  integration 
of  the  shear  force  or  moment  is  essential  in  the  first-order  forms. 

3.1.2  The  Bernoulli-Euler  Beam  Model 

The  Bernoulli-Euler  model  is  obtained  from  the  Timoshenko  model  with  two  additional 
assumptions: 
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1.  rotational  inertia  is  neglected,  pi  —*  0 

2.  shear  deformation  is  neglected,  —  <f>  — *  0 
Assumption  1.  reduces  the  second  of  equations  (54)  to 

■™(M -£("£)■  ,60) 

Equation  (60)  and  assumption  2.  are  now  used  to  reduce  the  first  equation  of  (54)  to 

.  <61> 


Note  that  equation  (60)  along  with  assumption  2  leads  to  the  following  expression  for  shear 
force 


(62) 


Although  (62)  is  commonly  used  in  conjunction  with  the  Bernoulli-Euler  model  (61),  it 
should  only  be  used  with  caution.  Equation  (61)  is  valid  only  in  the  limit  /  — »  0.  We  will 
return  to  this  point  below. 

The  boundary  conditions  (55)  and  (56)  reduce  to  for  a  =  0,  L 

displacement  or  shear  force 
•»(«,  «)  =  ».(«).  =  /„((),  (83) 

and 

displacement  or  moment 

^  =  *(<»,«)  =  *.(«).  =  »•.(<)  (64) 
Note  that  nonzero  shear  force  is  included  as  an  admissible  boundary  condition;  however, 
the  remarks  following  equation  (62)  apply. 

Equation  (61)  can  be  reduced  to  ‘first-order’  form  by  introducing  a  new  variable  7 (z,t) 

dr]  I  327 

dt 
dj 
dt 


Adz*' 

E  d2r] 


(65) 


(66) 


p  dz *' 

and  the  boundary  conditions  associated  with  (65)  are  for  a  =  0,  L: 

displacement  or  shear  force 

V(a,t)  =  »/a(t)>  =  7a(0»^  = 

and 

rotation  or  moment 

Observe  that  (65)  is  a  parabolic  system  of  the  type  (51).  Equations  (65)-(67)  can  be 
derived  directly  from  (61)  or  from  (57)  upon  invoking  assumptions  1  and  2.  We  should 
also  note  that  a  corresponding  expression  for  shear  force  becomes 


(67) 
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In  some  situations,  bending  deformation  may  be  negligible  with  respect  to  shear  deforma¬ 
tion,  that  is,  |^j  C  \dr)/dz\.  In  this  case,  the  first  equation  of  (54)  reduces  to 


pA 


d2r\ 
dt 2 


■  s  (-“£) 


(60) 


with  boundary  conditions  for  a  =  0,  L 

displacement  or 
q(a,t)  =  qa(t) 


shear  force 
kGA^  =  /.(*). 


(70) 


This  simple  model  is  primarily  useful  for  illustrative  purposes.  Again  by  introducing  the 
new  variable  u(z,t),  equation  (69)  can  be  replaced  by  two,  first-order  equations 


dr)  _  du 

dt  dz 
du  _  kG  dr) 

dt  p  dz 

subject  to  boundary  conditions  for  a  =  0,  L 

displacement  or  shear  force 

l(<M)  =  1..  «'(“>•)  *  *'•(0.  *’•  =  *$*. 


(71) 


(72) 


3.2  Beam  Models  with  Dissipation 

Various  dissipation  models  have  been  proposed  for  use  with  the  Timoshenko  and  Bernoulli- 
Euler  beam  models.  A  summary  of  the  most  frequently  cited  models  may  be  found  in 
Pich£  [53]  and  Wie  and  Bryson  [30].  Experimental  data  is  scant,  particularly  in  the 
high  frequency  range,  so  that  none  of  these  models  can  be  considered  definitive.  Recent 
experimental  work  by  Russell  has  demonstrated  that  the  usual  conjecture  of  ‘damping 
proportional  to  frequency’  for  internal  (material)  dissipation  is  not  complete  [82].  However, 
it  is  reasonably  representative  of  material  damping  for  certain  low  to  midrange  frequencies. 
Moreover,  some  of  the  more  popular  models  were  never  intended  for  vise  in  general  transient 
analysis  and  fail  to  yield  well-posed  dynamical  models  (Chen  and  Russell  [54]).  One 
approach  to  developing  dynamical  models  which  include  dissipation  is  to  augment  the 
variational  development  of  the  equations  of  motion  by  introducing  a  Rayleigh  dissipation 
function  (see  also  [82).  We  formulate  such  a  function  based  on  the  following  assumptions: 

1.  external  dissipation  forces  are  proportional  to  the  coordinate  velocities  (i.e.  17,^), 

2.  internal  dissipation  forces  are  proportional  to  strain  rates  (i.e.,  shear  strain  rate, 
If  =  dr)/dz  —  <f>,  and  compressive  strain  rate  d<f>/dz). 
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Cifp  rj  Kelvin- Voight  damping 

~crfcf  V  Chen-Russell  (or  ‘square  root’)  damping 
Ci*7  viscous  damping 


Table  1:  Dissipation  terms  and  standard  damping  mechanisms 
Thus,  we  define  the  dissipation  function  as 


(9ri\’  (  j  d  3„.  3*1’  (  3V \’l 

RM)  =  -Jo  +  «.(*)  +C3\ai(ar)'3r)  +  ’ 


where  Ci,C2  are  coefficients  of  external  damping  and  C3,e4  are  coefficients  of  internal  (ma¬ 
terial)  damping. 

The  modified  Timoshenko  equations  are  found  after  carrying  out  the  ususal  variational 
calculations  [27]  yielding 


.d7V  d  nA 

pAWt  =  r,  kGA 


thls-'I-’i'-IA-S) 

d<f>  (  d'r)  d<t>\  d3<t> 

Cj  dt  +  C3  \dzdt  dt  )  ^  C*dz7dt 


d3<t> 

'dz'dt 


To  obtain  the  Bernoulli-Euler  model  we  invoke  the  approximations  pi  -*  0  and  dr) / dz  — 
4>  — »  0  with  the  result 


52q  d4r) 

~  +  c»r?  +  EI^ZI 


pA"  +  “a?  ~ 


Table  1  summarizes  the  dissipation  terms  appearing  in  (75)  in  terms  of  standard  damping 
mechanisms  [51,54].  Individually,  these  dissipation  types  produce  well  known  eigenvalue 
patterns  for  the  spectrum  of  the  abstract  evolution  operator  A  in  (7).  In  particular, 
the  viscous  term  translates  the  spectrum  parallel  to  the  real  axis,  ‘square-root’  damping 
gives  the  wedge  pattern  characteristic  of  material  dissipation  and  Kelvin- Voight  damping 
provides  a  pattern  where  the  sequence  of  eigenvalues  have  an  accumulation  point  located 
on  the  finite  real  axis  [10].  See  for  example  Piche  [53],  Gibson  [10],  Chen  and  Russell  [54], 
Wie  and  Bryson  [30],  and  Balas  [9]. 

Equations  (73)-(74)  are  easily  put  in  the  first-order  form  appropriate  for  our  standard 
control  modeling  problem.  For  the  Timeshenko  model,  introduce  the  variables  i/(t,z), 
'littz)  and  write  the  equations  as 


dr) 

dv 

.  C3  d7r) 

dt 

dz 

pA  dz2 

du 

kG 

(On  A 

~di 

P 

C3  d<f> 


TT  ITT  —  c^. 


pA  dz 
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d<f>  _  di  (  A  c4  d2<t>  c3  dr\  (c2  +  c3) 

di  dz  I U  pi  dz2  pi  dz  pi  ’ 

dt  p  dz' 

For  the  Bernoulli-Euler  model  (75),  introduce  the  variable  7(2,*)  and  write  as 

dr? 
dt 
dj 
dt 

3.3  Frequency  Response  Calculations  for  DP  Systems 

In  this  section  we  will  be  concerned  with  the  computation  of  certain  irrational  transfer 
functions  and  a  resolvent  operator.  This  provides  a  complete  model  including  transient 
response  for  the  DP  system.  For  our  purposes,  the  resolvent  can  be  considered  as  an 
integral  operator  with  kernel  a  Green’s  Junction.  Using  transform  methods  we  compute 
explicit  formulae  for  the  abstract  objects  discussed  previously.  We  focus  on  hyperbolic 
and  parabolic  linear  (one  dimensional)  structural  models  for  distributed  elements.  Such 
models  can  be  used  for  elastic  dynamics  of  beams,  cables,  etc.  Finally,  the  computations 
are  extended  to  hybrid  system  models  consisting  of  interconnections  of  elastic  components 
with  rigid  bodies  and  other  lumped  parameter  models. 

3.3.1  Hyperbolic  Models 

Consider  a  class  of  elastic  structures  represented  by  hyperbolic  partial  differential  equations 
in  one  space  dimension  0  <  z  <  L,  (e.g.  arising  from  models  such  as  discussed  in  the 
previous  section): 

=  f  +  +  £„(,,,)  (78) 

subject  to  boundary  conditions 

Sjx(t,0)  -f  rii(t,L)  =  Df(t),  (79) 

and  initial  conditions 

z{Otz)  =  x°(x)  €  Mn(0,L).  (80) 

Here,  x  is  an  n-vector  valued  state  x  G  Hn(0,L),  v  G  Me(0 ,L)  is  an  ^-vector  valued  dis¬ 
tributed  disturbance,  /  is  m-vector  valued  boundary  interactions,  F,  H  are  real  n  x  n 
matrices  with  F  nonsingular  and  diagonalizable  [43],  and  Ei,  Tj  are  n  x  n  matrices.  Con¬ 
trollability  questions  for  systems  of  this  type  are  considered  in  Russell  (15].  We  remark 
that  (81)  is  a  concrete  example  of  the  transform  of  the  abstract  formula  (9).  After  taking 
Laplace  transforms  in  the  temporal  variable  t,  we  obtain 

=  jf  Gr(s,  z,  w)M  (a,  w)dw  +  Hbc(*,z)F(s) 


I  d2')  c4  d*r)  ci  d2rj 

A  A**  pA  dz*  +  pA  dz 2  ClT*’ 


Adz 2 
Ed2v 


(81) 


v*,\  ^  w,  >va  >.o  ar  w-  t«.v  sz  or  sag  r& 
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and  AT,  V,  F  are  the  Laplace  transforms  of  x,  t>,  /  respectively.  The  function  Gr(«,x,u;)  is 
the  Green’s  function  [50,52]  for  (78),  (70)  and  Hbc{*,*)  i*  a  transfer  function  from  bound¬ 
ary  interactions  to  state.  Since  in  most  cases  of  practical  interest  the  control  of  flexible 
structures  will  be  effected  by  actuators  whose  influence  functions  are  highly  localized,  we 
have  formulated  our  model  with  boundary  control  only. 

Comparison  with  (1)  clearly  shows  that  the  resolvent  for  the  operator  A  :  )/n(0,  L )  — ► 
^"(0,  L)  defined  by  (78)  and  (79)  is  the  integral  operator  /0L  Gr(s,z,w)  •  dxv. 

A  straightforward  calculation  leads  to  the  following  form  for  H&c 
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The  general  solution  is 

where  it  is  an  n-vector  to  be  determined.  Substituting  the  general  solution  into  the  bound¬ 
ary  conditions  yields 

Ei*  +  Txf-'W-Wk  =  Df[s) 

or 

k  =  [Ex  +  rieF"‘(w-*>I')-1£>/(s). 

By  substituting  k  into  the  general  solution  the  terms  in  (82)-(83)  can  be  recognized. 

Derivation  of  the  Green’s  function  follows  from  (84)-(85)  as  follows.  Note  that  inte¬ 
grating  (84)  over  a  small  interval  w  —  e  <  z  <  tv  +  t  gives 

rw+t  £  i-w+€ 

/  —Gr{s-,z,w)dz  =  /  F~1(sl  -  H)Gr{s\z,w)dz  +  /„ 

Jw—t  az  Jw—€ 

by  the  sifting  property  of  the  delta  function.  Then  by  letting  t  — ♦  0  we  see  that  the  second 
term  on  the  right  hand  side  vanishs  and  we  get 

Gr(s;tu+,u>)  -  Gr(s\w-,w)  =  (89) 

so  the  Green’s  function  is  discontinuous  at  z  =  w.  A  general  solution  to  (84)  has  the  form 
(86)  with 

GtLEft{s,z,w)  =  eF~,{‘I~H)*Kt 

GrRIGHT{s,Z,w)  = 

where  K(,Kr  are  n  x  n  matrices  to  be  determined  in  the  sequel.  Substituting  into  (89) 
gives 

Kr-Kt  =  c~F~  (90) 

and  into  (85)  gives 

E  iKt  +  VieF -'W-Wk,  =  0.  (91) 

Substitute  Kr  from  (90)  into  (91)  and  solve  for  Kt 

Kt  =  -[E, 

Now  substituting  Kt  into  (90)  and  solving  for  Kr  gives 

Kr  =  c-F-'l'I-H)*  -  (£,  +  rief-1W-H)i]-lrieF-‘(./-/r)(i-W) 

=  [/  -  [E]  + 

=  [e,  +  r,ef  [e,  +  r1eF",<w-">t  -  e'JP"1(*,-*)» 

=  [E,  +  r1er~,<*,-Jir>1)-IE, 

Substitution  then  into  the  general  solution  gives  the  forms  in  (87)-(88). 


SEI -TR-86-14 


29 


3.3.2  Parabolic  Models 

We  begin  with  the  canonical,  first-order  model  consisting  of  n  equations  as 

— ^ ^  +  g*(M)  +  S»(M)  (92) 

subject  to  2n  boundary  conditions 

Six(t,o)  +  +  rn(t,L)  +  =  P/(0,  (m) 

and  n  initial  conditions 

x(0,z)  =  x°(z)  e  Un{0,L).  (94) 

Let  E  =  [Ei,E2]  and  T  =  [ri,T2] — each  2n  x  2n  real  matrics. 

Following  a  procedure  as  above  the  Green’s  function  and  boundary  transfer  function 
can  be  computed  as  follows.  Let 


A(fi)  =  -G~'{H-sIn)  -G~lF  ’ 

9{s,z)  =  eA<'>4, 

M(s,z)  =  (Jn,0]$(s,s)(£  +  r$(s,.L)]-1. 


Then  the  boundary  transfer  function  is 


HBc{s,z)  =  M(s,z)D, 


and  the  Green’s  function  is  given  by 


„  /  %  _  (  GtLeft{s,z,w),  forO  <  z  <  w 

1  GrRtGHT{s,z,w),  lot  w  <  z  <  L 


GrlEFT  =  -Af(j,a)r#(«,L -w)  °  , 

+  n 

Gr  RIGHT  =  Af(s,x)E$(s, -Itf)  J 

*n 


3.4  Modeling  of  Hybrid  Systems 

In  most  applications,  models  for  the  dynamics  of  flexible  structures  will  involve  interaction 
between  various  elastic  and  rigid  elements.  In  the  particular  case  of  flexible  structures 
associated  with  large  space  structures,  the  potential  topological  configurations  can  be 
quite  complex.  Various  elements  such  as  beams,  truss  structures,  cables,  membranes,  etc., 
may  have  dominant  distributed  parameter  effects.  Typically  a  central  body  or  bodies 
represent  large  concentrations  of  mass  with  respect  to  the  overall  low  mass  density  of  the 
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flexible  structure.  These  are  most  effectively  represented  by  lumped  parameter  models 
of  their  rigid  body  dynamics.  Additionally,  various  attitude  control  actuators  can  add 
concentrated  inertia  elements  which  can  be  effectively  modeled  as  lumped  systems.  Thus, 
carefully  chosen  linear,  hybrid  models  can  provide  an  effective  tool  for  analysis  of  dynamics 
of  vibrations  and  their  effect  on  small  angle  motions  for  complex  space  platforms.  In  this 
section,  we  consider  the  structures  and  computations  of  certain  resulting  transfer  functions 
and  the  resolvent  operator  for  the  composite  system  along  the  lines  of  Section  1. 

The  concept  of  a  mechanical  impedance  (terminology  borrowed  from  electrical  network 
theory)  has  been  used  in  structural  dynamic  modeling  for  many  years  [51].  The  dynamic 
stiffness  method  (application  to  space  structure  modeling  is  reviewed  in  Piche  [53])  uses 
this  notion  to  compute  effective  transfer  function  models  for  interconnected  structures  [55]. 
Our  approach  here  will  follow  along  similar  lines  except  that  we  will  focus  on  computing  the 
resolvent  operator  for  a  hybrid  structure  by  direct  manipulation  of  its  kernel;  viz,  a  Green’s 
function  [52].  A  hybrid,  state  space  model  is  constructed  in  Burns  and  Cliff  [12]  (where 
considerations  are  given  for  approximation  and  computation  in  the  hybrid  state  space). 
We  will  consider  a  hybrid  state  space  as  consisting  of  a  direct  sum  of  spaces  X  =  X t  ©  Xd 
where  Xd  =  is  the  distributed  part  constructed  on  an  appropriate  Hilbert  space  of  JVj- 
vector  valued  functions  with  the  “energy”  inner  product  of  (6)  for  a  distributed  parameter 
system  (DPS)  written  in  abstract  form  (7)  and  Xe  =  fRNl,  a  finite  dimensional  state  space 
of  the  lumped  parameter  system  (LPS). 

For  control  of  hybrid  structures,  we  restrict  attention  to  DPS  modeled  as  either  the 
hyperbolic  or  parabolic  (or  mixed)  cases  which,  as  we  have  seen,  can  be  expressed  in  the 
frequency  domain  in  the  form 


rL 

_Xd(s,z)  =  /  Gr(s,z,xv)M(s,w)dw  +  HBc[s,q)Fd{s) 
Jo 


(100) 


N 

> 

> 

> 

\ 


3 


8 


*  • 
S. 


where 

M(s,w)  =  x°(u;)  —  EV  (s,  w).  (101) 

Clearly,  (100)-(l0l)  can  represent  a  disjoint  collection  of  distributed  elements  such  as 
beams,  cables,  etc.  (Conceptually,  a  version  of  (100)  can  also  be  written  for  higher  dimen¬ 
sional  spatial  domains,  but  we  feel  for  the  current  presentation  that  the  required  complexity 
of  notation  can  mask  the  simplicity  of  the  underlying  concepts — see  Butkovskiy  for  details 
[52].) 

All  LPS  component  models  are  combined  into  a  LPS  state  space  model  as 

ie(t)  =  Aexe(t)  +  Btft{t),  x\  =  i/(0)  (102) 

with  x/  €  RNl  =  Xt  a  finite  dimensional  real  space.  By  taking  Laplace  transforms  in  (102), 
we  write  (analogous  to  (100)) 


I 


4 


R 


Xt(s)  =  Ri(s)if  +  Ht(s)Fe(s), 


(103) 


where  i?/(s)  =  [sJ/yr,  —  A/]-1  is  the  resolvent  for  the  (matrix)  operator  At  and  Bt(s)  = 
Rt{s)B. 
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The  hybrid  state  space  X  =  Xt  ©  X  d  consists  of  elements 


*■■>  ■  ( ) 


1.104) 


which  are  N  =  JVd  +  AV valued  functions  of  z  6  [0,  L],  t  >  0.  Finally,  the  interconnection 
of  component  systems  is  resolved  through  a  topological  constraint  relation  consisting  of 
m  =  rrid  +  m*  linear  equations; 


f(t)  +  7|Zrf(f»0)  +  T2Xd(t,L)  -f-  T$xt(t)  —  Ku(t) 


(105) 


where  u(t)  is  a  fc-vector  of  control  inputs  to  the  hybrid  system,  T(,T2  are  m  x  Nd,  T3  is 
m  x  Nt ,  and  K  is  mxk  real  matrices.  The  hybrid  modeling  problem  is  to  find  an  equation 
of  the  form  (100)  by  solving  (100),  (103)-(105)  simultaneously  for  the  hybrid  state  x[t,z). 
We  provide  the  resulting  model  in  the  following  form: 


AT(s,z)  =  /  Gr(s,z,tv)M(8,w)dw  +  R{s,z)x°t  +  Hbc{«,z)U{s), 

Jo 

where  M(s,u>)  is  given  in  (101).  The  resolvent  operator  for  the  hybrid  system  is 

JZ(s;A)  =  R[s,z),  f  Gr(s,z,w)  •  dw 

Jo 


(106) 


(107) 


where  R(s;  A)  :  X  -*  Z>(A)  C  X ,  Gr  is  IV  x  and  R  is  N  x  JV*  are  matrix  valued  functions 
which  can  be  computed  explicitly  as  follows: 


Gr(s,z,ti;)  = 


~#/(s)Qi(s) 

Gr{s,z,w)  -  Hbc(s,z)Q2(s ) 


P{s,w) 


where 


Q(s)  =  1  /„  +  Q(s)i-  =  [  ' , 

Q{«)  =  \T3Ht{s),  T\Hbc{s,0)  +  T2Hbc{s,L)\  , 

P{s,w)  =  TiGr(s,0,io)  +  T2Gr(s,L,w). 


(108) 


(109) 


(110) 

(111) 

(112) 


Finally,  the  N  x  it  transfer  function  matrix  from  boundary  control  to  hybrid  state  is 


_  f  «<(») 


0  Hbc{s,z) 


Q(s)K. 


(113) 


The  derivation  of  (106)-(112)  is  straightforward  and  proceeds  as  follows.  Substitute 
(100),  (102)  into  (105)  and  solve  for  the  interconnecting  force  F(s).  This  identifies  the 
terms  Q(»),P(stw)  above.  Now  substitute  the  appropriate  components  of  F(s)  into  (100), 
(102)  and  use  the  hybrid  state  model  (104). 

In  Section  5  of  this  report  we  consider  several  examples  illustrating  these  calculations. 
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4  Homogenization  of  regular  structures 

In  this  section  we  consider  the  problem  of  modeling  and  control  of  a  class  of  lattice  struc¬ 
tures,  e.g.,  trusses.  Their  large  size  and  repetitive  infrastructure  require  special  techniques 
for  structural  analysis  to  cope  with  the  large  number  of  degrees  of  freedom.  Approxi¬ 
mations  of  such  systems  by  continua  provide  a  simple  means  for  comparing  structural 
characteristics  of  lattices  with  different  configurations,  and  they  are  effective  in  repre¬ 
senting  macroscopic  vibrational  modes  and  structural  response  due  to  temperature  and 
load  inputs.  Our  approach  to  the  construction  of  such  models  is  based  on  a  technique 
for  asymptotic  analysis  called  homogenization.  It  has  been  widely  used  in  mathematical 
physics  for  the  treatment  of  composite  systems  like  porous  media  for  which  one  wishes 
to  have  an  effective  approximating  system  with  parameters  which  are  constant  across  the 
structure.1 

Before  developing  the  general  features  of  the  method  and  applying  it  to  the  treatment 
of  lattice  structures,  we  shall  make  a  few  remarks  on  other  work  on  continuum  models 
which  has  appeared  in  the  recent  structural  mechanics  literature. 

Noor,  et.  al.  [l]  use  an  energy  method  to  derive  a  continuum  approximation  for  trusses 
with  triangular  cross  sections  in  which  the  modal  displacements  of  the  truss  are  related  to  a 
linearly  varying  displacement  field  for  an  equivalent  bar.  Plates  with  a  lattice  infrastructure 
are  also  treated.  In  Dean  and  Tauber  [64]  and  Renton  [74]  exact  analytical  expressions 
for  the  solutions  of  trusses  under  load  were  derived  using  finite  difference  calculus.  By 
expressing  the  difference  operators  in  terms  of  Taylor’s  series  Renton  [77]  was  able  to  derive 
continuum  approximations  to  the  finite  difference  equations  resulting  in  expressions  for 
equivalent  plate  stiffnesses,  for  example.  In  a  recent  paper  Renton  [77]  used  this  approach 
to  give  equivalent  beam  properties  for  trusses,  which  complements  the  earlier  work  of  Noor, 
Anderson  and  Greene  [3],  and  Nayfeh  and  Hefzy  [2].  (See  also  (Anderson  [3]).) 

In  most  cases  a  continuum  model  is  associated  with  the  original  (lattice)  structure 
by  averaging  the  parameters  of  the  lattice  over  some  natural  volume  (e.g.,  of  a  “cell”  of 
the  structure)  and  identifying  the  averaged  parameter  value  (mass  density,  stress  tensor, 
etc.)  with  the  corresponding  distributed  parameter  in  the  continuum  model.  A  specific 
form  for  the  continuum  model  is  postulated  at  the  outset  of  the  analysis;  e.g.,  a  truss  with 
lattice  structure  will  be  approximated  by  a  beam,  with  the  beam  dynamical  representation 
assumed  in  advance.  While  this  approach  has  an  appealing  directness  and  simplicity,  it 
has  some  problems. 

First,  it  is  very  easy  to  construct  an  example  in  which  the  “approximate  model”  ob¬ 
tained  by  averaging  the  parameters  over  a  cell  is  not  a  correct  approximation  to  the  system 
behavior.  This  is  done  in  subsection  4.1. 2  Second,  one  cannot  use  this  procedure  to  obtain 
“corrections”  to  the  approximation  based  on  higher  order  terms  in  an  expansion,  which 
may  sometimes  be  done  in  an  asymptotic  analysis.  These  terms  can  be  used  to  describe  the 
microscopic  behavior  (e.g.,  local  stresses)  in  the  structure.  Third,  the  averaging  method 


'See,  for  example,  the  paper*  of  Lanen  [89],  Keller  [68],  and  the  report*  of  Babueka  [57]  for  application* 
and  discussion*  of  design  technique*. 

3See  the  numerical  experiment*  in  (Bourgat  [62]). 


(averaging  the  parameters  over  space)  does  not  apply  in  a  straightforward  way  to  Bystems 
with  a  random  structure,  since  the  appropriate  averaging  procedure  may  not  be  obvious.3 
Fourth,  the  method  cannot  be  naturally  imbedded  in  an  optimization  procedure;  and  con¬ 
trols  and  state  estimates  based  on  the  averaged  model  may  not  be  accurate  reflections  of 
controls  and  state  estimates  derived  in  the  course  of  a  unified  optimization  -  averaging 
procedure.  In  particular,  the  method  does  not  provide  a  systematic  way  of  estimating 
the  degree  of  suboptimality  of  controls  and  state  estimates  computed  from  the  idealized 
model. 

In  this  work  we  use  a  totally  different  technique  called  homogenization  from  the  math¬ 
ematical  theory  of  asymptotic  analysis  to  approximate  the  dynamics  of  structures  with  a 
repeating  cellular  structure.  Homogenization  produces  the  distributed  model  as  a  conse¬ 
quence  of  an  asymptotic  analysis  carried  out  on  a  rescaled  version  of  the  physical  system 
model. 

Unlike  the  averaging  method,  homogenization  can  be  used  in  combination  with  opti¬ 
mization  procedures;  and  it  can  yield  systematic  estimates  for  the  degree  of  suboptimality 
of  controls  and  estimators  derived  from  idealized  models.  While  our  results  are  stated 
in  terms  of  simple  structures,  they  demonstrate  the  feasibility  of  the  method;  and  they 
suggest  its  potential  in  the  analysis  of  structures  of  realistic  complexity. 

In  subsection  4.1  we  give  an  example  derived  from  (Bensoussan,  Lions,  and  Papanico¬ 
laou  [60])  illustrating  some  of  the  subtleties  of  homogenization,  particularly  in  the  context 
of  control  problems.  In  subsection  4.2  we  derive  a  homogenized  representation  for  the  dy¬ 
namics  of  a  lattice  structure  undergoing  transverse  deflections.  We  show  that  the  behavior 
of  the  lattice  is  well  approximated  by  the  Timenshenko  beam  equation;  and  we  show  that 
this  equation  arises  naturally  as  the  limit  of  the  lattice  dynamics  when  the  density  of  the 
lattice  structure  goes  to  infinity  in  a  well  defined  way.  The  problem  of  vibration  control  of 
a  lattice  is  posed  and  discussed  in  subsection  4.3.  In  subsection  4.4  we  derive  a  diffusion 
approximation  for  the  thermal  conductivity  of  a  one-dimensional  lattice  structure.  This 
property  is  useful  in  analyzing  new  materials  for  large  space  structures. 

Acknowledgements:  We  are  grateful  to  Professor  George  Papanicolaou  for  bringing 
Kunnemann’s  paper  to  our  attention  and  to  Drs.  A.  Amos  and  R.  Lindberg  for  their 
comments  on  an  earlier  version  of  this  work. 

4.1  A  one-dimensional  example 

From  (Bensoussan,  LionB,  and  Papanicolaou  [60])  we  have  the  following  example: 

~  =  /(x)»x  e  (xo,*i) 

u‘(z0)  =  0  =  u‘(zi) 

’Homogenization  methods  do  apply  to  systems  with  a  randomly  heterogeneous  structure,  see  (Papan¬ 
icolaou  and  Varadhan  [71])  and  (Kunnemann  [78]).  We  shall  treat  such  systems  in  a  subsequent  report. 


(114) 
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where  o‘(x)  =f  a(x/e),  and  a(y)  is  periodic  in  y  with  period  Yo>  «(v)  >  a  >  0.  It  is  simple 
to  show  that 

ll-'lli.  =  jT‘  KWI-  +  <  C  (115) 

and  so,  u*  — ►  u  weakly  in  the  Hilbert  space  ff1.4  Moreover, 

a*  -♦  M{a)  =f  f  a(y)dy  (116) 

I o  •'0 

and  it  is  natural  to  suppose  that  u*  — ►  u  with  the  limit  defined  by 

~  =  f(x)’x  €  (*o,*i)  (117) 

u(x0)  =  tt(x,) 

This  is  untrue  in  general  (Bensoussan,  Lions,  and  Papanicolaou  [60],  pp.  8-10).  The 
correct  limit  is  given  by 

“  ^l6^u(z)l  =  /(*)» x  e  (xo,  *i)  (118) 

u(x0)  =  u(x,) 

with 

a^[M(i))-’  (119) 

In  general,  M(a)  >  a;  and  so,  the  error  in  identifying  the  limit,  (117)  versus  (118),  is 
fundamental. 

The  system  (117)  corresponds  to  averaging  the  parameter  a'(x)  over  a  natural  cell;  a 
procedure  similar  to  that  used  in  the  past  to  define  continuum  models  for  lattice  structures. 
As  (118)  shows,  the  actual  averaging  process  can  be  more  subtle  than  one  might  expect, 
even  for  simple  problems. 


4.1.1  Homogenisation  of  the  example 

To  see  how  (118)  arises,  we  can  use  the  method  of  multiple  scales  which  applies  to  a  variety 
of  perturbation  problems.  Suppose 

u€(x)  =  u*(x,  *)  =  uq(x,  *)  +  etii(x,  *)  +  . . .  (120) 


that  is,  we  suppose  that  u*  depends  on  the  “slow"  scale  x  and  the  “fast”  scale  y  =f  x/c; 
and  we  adopt  an  ansatx  which  reflects  this  dependence.  Using  the  identity 


S'*- 7» 


du  ldu  x 

T  +  “3“*^  =  “ 
ox  t  ay  e 


•Her.  «'  t'  {.  6  £>(«.,.,) :  ||.||H.  <  oo) 


(121) 
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then  (114)  may  be  rewritten  as 

,  d  1  d  Vf  ,  w  d  Id,.  ,, 

~{ai  +  7^'{“(v)(ai  +  Ia^“°  +  tu‘ +  -1>  =  > 

Simplifying  and  equating  coefficients  of  like  powers  of  t,  we  find  first  that 

1  d  ,  .  ,  d  , 

"  m  °- 

The  assumptions  on  a(y)  imply 

«o(*,y)  =  u0(x) 

i.e.,  no  y-dependence.  The  coefficients  of  «_1  satisfy 


If  we  look  for  «i  in  the  form 


d  ...  d  .  dadua 


«i(x»y)  =  -x(y)-^  +  Mi), 


then  the  corrector  x(y)  must  satisfy 


and  be  periodic.  That  is, 


-  — [a(y)— y(»)l  = 

dy  ^  dy  dy 


o(y)^  =o(y)  +  e 


(122) 


(123) 


(124) 


(125) 


(126) 


(127) 


(128) 


(129) 


which  has  a  periodic  solution  (unique  up  to  an  additive  constant  in  y)  if  and  only  if 


which  implies 


1  rr°  c 

n/o  '1  +  ^  =  ° 

«--Wr)|-®» 


(130) 


(131) 


We  obtain  an  equation  for  u0(x)  from  the  solvability  condition  for  u2(x,y).  Equating  the 
coefficients  of  e°  in  the  expansion,  we  have 


(132) 


.  «  L  »  il%_|lt)  *'«■  |k, 


,  4U  «!,  ‘“Jt-  ->«■-  ■■  .>  -.-« 
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This  has  a  solution  uj(x,y),  periodic  in  y  if  and  only  if 


{ y0J0  °(afr) +  ^[a(v)x(v)l  -  a(y)^x(y)]^y} 


dP\io{x) 
dx 2 


(133) 


+/(i)  =  0 


where  we  have  used  (127).  The  integral  of  the  second  term  is  zero,  since  it  is  the  integral  of 
the  derivative  of  a  periodic  function  over  one  period.  Using  (129)  and  (131),  (133)  reduces 
to 

-4^  +  /(x)=0  .  (134) 


(plus  the  boundary  conditions)  which  is  (118)  (119). 


4.1.2  Control  and  homogenisation  of  the  one  dimensional  system 


One  of  the  simplest  stochastic  control  problems  associated  with  the  preceding  system  is 
defined  by  the  Hamilton  -  Jacobi  -  Bellman  equation 


.X'.dPu*  1,-i.dti* 

-“Vis  ~  d* 


(135) 


l€0,u'(z)  =  Oonr  ='  SO 

where  O  is  an  open  interval  in  R,  and  each  function  a(y),b(y),  and  y(x,y)  is  periodic  in 
y  with  period  y0.  We  assume  that  a(y)  >  a  >  0  and  that  c  >  0,  and  that  the  controls  v 
take  values  in  R. 

This  Bellman  equation  corresponds  to  the  stochastic  control  problem 


**'(*)  =  l") 


•/>(.)]  =  Em{P  £<*'*'>*]*} 

Jo  e 


(136) 


dx‘(t)  =  o(x‘,  )dw(t )  +  -fe(x',  —)dt  -f-  G(xt,  —,v)dt 


x‘(0)  =  xeo,i>o. 

with  <r2(x,y)  =  a(y),  6(x,y)  =  6(y),  G(x,y,v)  =  ff(i,y)v,  /(x,y,v)  =  §v2,  and  c(x,y, v)  = 
c,  a  constant  in  (135).  Each  function  in  (136)  is  assumed  to  be  periodic  in  y  with  period 
one.  We  are  interested  in  the  behavior  of  the  optimal  cost  and  control  law  for  (135)  in 
the  limit  as  c  — ♦  0.  The  stochastic  control  problem  (136)  was  treated  in  (Bensoussan, 
Boccardo,  and  Murat  [59]);  the  analysis  here  uses  different  arguments  which  emphasize 
the  computational  aspects  of  the  system. 


>.*•>  >?m  9.  o;<  <“vs  ^  % 
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(139) 


(140) 


Evaluating  the  infimum  in  (135),  we  have  the  nonlinear  system 

“(f  K.  +  i*(*K  -  «■',  - -t  )(«*,)'  =  o  (137) 

x  e  0,uf(x)|r  =  0. 

The  analysis  of  the  control  problem  involves  homogenization  of  this  system. 

Let 

M  -  a(y)dvv  +  b(y)dv  (138) 

with  its  formal  adjoint  defined  by 

A\  =  av[a(y)^-]  -  3w((6(y)  -  a„(y))-].  *  (139) 

The  problem 

A\m  =  0,y  -»  m(y)  periodic  (140) 

m  >  0,  L  m(y)dy  =  1 

has  a  unique  solution  m(-)  onK  =  5°,  the  unit  circle,  with 

0  <  m  <  m(y)  <  m  <  oo.  (141) 

So  m(-)  is  a  density  on  Y.  We  assume  that  fc(-)  is  centered 

JY  ™{y)Hy)dy  =  o.  (142) 

As  a  consequence  the  system 

Aix{y)  =  b{y)  (143) 

y  -*  x(y)  periodic  ,J  x{y)dy  =  0 

has  a  well  defined  solution.  x(*)  is  the  corrector  associated  with  the  problem. 

As  before  we  set  y  =  x/e  and  look  for  u‘  in  the  form 

u*(x)  =  u‘(x,y)  =  u0(x,y)  +  €Ui(x,y)  +  . . . ,  (144) 

and  we  use 

d,<f>{x,  y)  =  <£x(x,  y)  +  - 4>v{x ,  y),  y  =  x/e  (145) 

d*x<t>{x,y)  —  <£,*(x,y)  +  2c<fl  «v(z»y)  +  &w(x>  y)- 
Substituting  in  (137),  we  have 

a(y)[uo«  +  2eu0«v  +  +  a(y)[fu,«  +  2uj«„  + 

+a(y)(eJu2„  +  2(u2ty  +  u2|fV] 

+  j6(y)[u0«  +  ^uoy)  +  j6(y)[cuu  +  uI|f]  (146) 


(141) 


(142) 


(143) 


(144) 


(145) 
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-t-j6(y)[f2u2s  +  £u2y]  -  c[u0  +  eui  +  £2tt2] 

-4y2(*>y)  [(«<>*  +  culs  +  c2u2,)  +  - (uqj,  +  fu,„  +  £2u2y)]2  =  0(£2) 
i  € 

The  last  term  is  ^ 

-  rSJ(a:»y)[(-2«o»  +  2eu0zuoy  +  2u0,uiy  +  u^) 

+c(2uo*Ux*  +  2ulxui„  +  2uiyu2y  +  2uo*u2y)]  +  0(£2) 

Equating  coefficients  of  like  powers  of  e,  we  obtain 

(c“2)a(y)uovw  +  *>(y)«ov  -  ^sJ(*>y)«o»  =  0  ' 

(£-1)a(y)ulyv  +  b{y)uXy  +  2a(y)u0„ 

+6(l/)uo*  ~  yJ(^>y)u o*«ov  =  0 
(t°)a(y)u2yy  +  fc(y)u2y  +  2  a(y)uley  +  b[y)uu 

+a(y)u0„  -  cu0  -  ^ff2(x,y)u^  -  y2(x,y)u0x«i*  =  0. 

Choosing  u0(z,y)  =  u0(z),  which  must  be  justified,  satisfies  (148).  We  can  then  solve  (149) 
by  choosing 

ui(*.y)  =  ~x(y)“o*(x)  +  ux(z).  (151) 

Equation  (150)  has  a  solution  for  u2(x,y)  if 

J  m(y){-2a(y)xwu0«  -  &(y)Xy“o*  +  a(y)«o«  (152) 

-c« o  -  ^y2(z,y)(l  -  2xMx}dy  =  0. 

This  gives  an  equation  for  uq(z) 


(147) 

(148) 

(149) 

(150) 


4«o..  -  cuo  -  ^ru;,  =  0 


where 


9  =f  JY  m(y){fl(y)[i  -  2xy(y)l  -  x(y)&(y)}dy 
r  ^  Iy  “  2Xy(y)]<*y- 

Remark.  From  the  definition  of  Ax  and  the  corrector  x(y)  we  have 

Jy  m(y)6(y)x(y)dy  =  Jy[a{y)xn  +  &(y)x*]x(yMy)dy 

=  jy  X(y)5vv[a(y)x(y)m(y)]dy  -  x[v)dv\b{y)x{y)m(y)}dy 


(155) 
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Also,  using  (143), 

fy  m(y)b{y)x{y)dy  =  x(y)a(vMy)Xw,  {v)dy  (156) 

-  J  x{v)Hv)m(v)x*dV  +  2Jy  x(y)XiA[o(s/)m(y)]dy 

Adding  these  two  expressions,  we  have 

2  fYm{v)b{v)x{v)dV 

=  2  JYx{y)*{y)m[y)xwdy  + 2  fYx{y)xv[am}vdy  (157) 

=  -2  JY  dy\xam]xvdy  +  2  j  x{y)xvl*rn]vdy  -  2  J^Xva{y)m{y)xvdy 
Thus,  q  may  be  rewritten  as 

«  =  ^  m(y){a(y)l1  -2Xv  +  x\\)dy  (158) 

=  JYm(y)a{y)ll  ~  xv?dy 

and  clearly  q  >  0. 

The  term  q  in  (154)  summarizes  the  effects  of  the  averaging  process  on  the  uncontrolled 
system.  The  homogenization  process  interacts  with  the  control  system  through  the  term 
T,  whose  form  would  be  difficult  to  “guess”  from  simple  averaging  procedures. 

4.2  Continuum  Model  for  a  Simple  Structural  Mechanical  Sys¬ 
tem 

4.2.1  Problem  definition 

Consider  the  truss  shown  in  Figure  1  (undergoing  an  exaggerated  deformation).  We  shall 
assume  that  the  truss  has  a  regular  (e.g.,  triangular)  cross-section  and  no  “interlacing” 
supports.  We  assume  that  the  displacements  of  the  system  are  “small”  in  the  sense  that  no 
components  in  the  system  buckle.  We  are  interested  in  describing  the  dynamical  behavior 
of  the  system  when  the  number  of  cells  (a  unit  between  two  (triangular)  cross  sections)  is 
large;  that  is,  in  the  limit  as 

c  =f  l/L  ->  0.  (159) 

We  shall  make  several  assumptions  to  simplify  the  analysis.  First,  we  shall  assume  that 
the  triangular  sections  are  essentially  rigid,  and  that  all  mobility  of  the  system  derives 
from  the  flexibility  of  the  members  connecting  the  triangular  components.  Second,  we 
shall  ignore  damping  and  frictional  effects  in  the  system.  Third,  we  shall  confine  attention 
to  small  transverse  displacements  rf(t,x)  and  small  in  plane  rotations  ^(f,x)  as  indicated 
in  Figure  1,  ignoring  longitudinal  and  out  of  plane  motions  and  torsional  twisting.  Fourth, 


j* 
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i) 


Figure  1:  Deformed  truss  with  regular  cross-section. 

we  shall  assume  that  the  mass  of  the  triangular  cross  members  dominates  the  mass  of  the 
interconnecting  links. 

Systems  of  this  type  have  been  considered  in  several  papers  including  [l],  [2],  {3}, 
and  [77].  In  those  papers  a  continuum  beam  model  was  hypothesized  and  effective  values 
for  the  continuum  system  parameters  were  computed  by  averaging  the  associated  param¬ 
eters  of  the  discrete  system.  Our  approach  to  the  problem  is  based  on  homogenization  - 
asymptotic  analysis  and  is  quite  different. 

The  assumptions  simplify  the  problem  substantially,  by  suppressing  the  geometric 
structure  of  the  truss.  We  can  retain  thi6  structure  by  writing  dynamical  equations  for  the 
nodal  displacements  of  the  truss  members.  For  triangular  cross  sections  nine  parameters 
describe  the  displacements  of  each  sectional  element.  The  analysis  which  follows  may  be 
carried  over  to  this  case,  but  the  algebraic  complexity  prevents  a  clear  presentation  of  the 
main  ideas.  As  suggested  in  (Noor  et  al.  [l])  one  should  use  a  symbolic  manipulation 
program  like  MACSYMA  to  carry  out  the  complete  details  of  the  calculations.  We  shall 
take  up  this  problem  on  another  occasion;  for  now  we  shall  treat  the  highly  simplified 
problem  which,  as  we  shall  see,  leads  to  the  Timoshenko  beam. 

We  shall  begin  by  reformulating  the  system  in  terms  of  a  discrete  element  model  as 
suggested  in  (Crandal  et  al.  [27]);  see  Figure  2.  In  this  model  we  follow  the  displacement 
and  rotation  <f>i{t)  of  the  itk  mass  M.  The  bending  springs  (k\)  tend  to  keep  the 
system  straight  by  keeping  the  masses  parallel  and  the  shearing  springs  (k\)  tend  to  keep 
the  masses  perpendicular  to  the  connecting  links.  We  assume  small  displacements  and 
rotations  so  the  approximations 

sin &(f)  =*&(*)  (160) 

tan-1  MO/*]  =  »fc(0/* 

are  valid. 


.•>  ••  *»  . 
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Normalizing  t  =  1,  we  associate  a  position  x  E  [—1/2, 1/2]  with  each  mass;  and  we  intro¬ 
duce  the  notation 

ri(t,Xi)  =  rii(t),4>(ttXi)  =  <£,(*).  (166) 

Having  normalized  t  —  1,  we  have  c  =  t  and  x<+ \  —  x,-  + 1  —  x,  +  c.  Let  Z  =  {x,}  be  the 
set  of  all  points  in  the  system.  In  this  notation 

(Ve+rj)(t,x)  =  ^[f?(*,x  +  c)  -  rj(t,x)]  (167) 

(vt~»?)(*ia:)  =  ^  [»?(<.*)  “  *?(*.*  -  e)],x  e  Z 

and  the  system  is 

=  K,(zi){V'*T,'(t,xi)  -  *'(1,1,)} 

+  rV'+{Kl(ij)V-*'(t,i,))  (168) 

<P’>^’Ii)  =  —  V-{.K.(x,)[V'V(l,ii)  -  *‘(1, *()]},*  e  z 

The  scaling  of  (168)  may  be  interpreted  in  the  following  way:  Formally,  at  least,  the 
right  sides  of  both  terms  in  (168)  are  0(t_a).  This  implies  that  the  time  variations  are 
taking  place  in  the  “fast  time  scale”  r  =  t/e.  Also,  the  spatial  variations  are  taking  place 
in  the  “microscopic  scale”  x  which  varies  in  c-increments  (e.g.,  x,+i  =  x,  +  c).  Introducing 
the  macroscopic  scale  z  =  ex,  and  the  slow  time  scale  a  —  cr,  we  may  rescale  (168)  and 
observe  its  dynamical  evolution  on  the  large  space-time  scale  on  which  macroscopic  events 
(e.g.,  “distributed  phenomena”)  take  place. 

Rewritten  in  this  spatial  scale,  the  system  becomes 

-  **(!.*<)} 

+  i«'+{rjr,(V)<-*'(t,*')}  (169) 

=  ^-(a:.(V)[«,+»((,V)  -  <*•(!,  *•))}  (170) 

where 

S'*  =  £VC±  =  0(1)  in  c.  (171) 

The  essential  mathematical  problem  is  to  analyze  the  solutions  of  (169)(170)  in 
the  limit  as  e  — *  0. 

4.2.2  Mathematical  analysis 

To  proceed,  we  shall  generalize  the  problem  (169)  (170)  slightly  by  allowing  Kt  and 
to  depend  on  z  as  well  as  zje.  This  permits  the  restoring  forces  in  the  model  system  to 
depend  on  the  large  scale  shape  of  the  structure  as  well  as  on  local  deformations.  We 
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(172) 


use  the  method  of  multiple  scales;  that  is,  we  introduce  y  =  x/e  and  look  for  solutions  of 
(169)  (170)  in  the  form 

»7*(0  =  v*{t,  *.y)»4'(0  =  MM.y)* 

and  we  have 

if,  =  K,(z,y),Kt  =  Kt(z,y),y  =  ^ 

On  smooth  functions  ^(z,  zf)  the  operators  £f±  satisfy 

(6e+V»)(2,y)  =  rl>(z  +  f,y  +  1)  -  tp{z,y) 

=  V»(2,y  + 1)  -  V>(*.y)  +  V»(*  +  e,y  + 1)  -  ^(*,y  + 1) 


(174) 


=  (5+^)(z,y)  +  «^j(*,y  +  1)  +  +  *)  +  0(f3) 


(6*"0)(z,y)  =  V>(*.y)  -  V>(*  -  £,y  -  i) 

=  ^(*,y)  -  V>(*,y  -  l)  +  V>(*,y  - 1)  -  V>(*  -  £.y  -  i) 


=  (5  tl>)[z,y)  -  e§^(*.v  +  l)  +  |«*§^f(*.v  ~  +  °(f3) 


(175) 

(176) 


We  assume  that  <f> *  and  rjf  may  be  represented  by 

f(<.«,y)  =  MM)  +  «MM,y)  +  ... 

»?€(t,z,y)  =  rj0(t,z)  +  e»?i(M>y)  +  ... 
and  substituting  (176)  in  (169)(170)  and  using  (173)  (174) (175),  we  arrive  at  a  sequence 


of  equations  for  (<Ao,»/o)>  (M*h)>  •  •  .  by  equating  the  coefficients  of  like  powers  of  e. 

r-2  ,o 


Starting  with  e~J,  t-1,  t°, . . .,  we  have 


-S+[rtf(z,y)S-,Mt,2)]=0 


(177) 


which  is  trivially  true  from  (175)  (176).  The  same  term  involving  qo(M)  from  (170)  is 
trivially  satisfied  by  the  assumption  (176).  Continuing 


7l5+{r-M*»y)5~MM,y)} 


(178) 


+K.{z, y){S+f?0(f,z)  -  MM)}]  =  0 
which  may  be  solved  by  using  the  corrector  x+(z,  y)  and  taking 

MM.y)  =  x#(*,y)MM) 

with 

S+{r.Kt(z,y)S_x*(*,y)}  =  K,(z,y) 

If  we  regard  z  as  a  parameter  in  (180),  then  there  exists  a  solution  x#>  unique  up  to  an 
additive  constant,  if  A*(z,  •),  Jf,(z,  •)  are  periodic  in  y,  if  there  exist  constants  A  and  B  so 


(179) 

(180) 


that 


0  <  A  <  Ki(z,y)  <  B  <  oo 


(181) 


. \ 
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and  if  the  average  of  K,(z ,  •)  is  zero 

1  rL/* 

t  L,7K-i’'V)dy  =  0 

Let  us  assume  that  (181)  (182)  hold,  and 

0  <  A  <  Kt(z,y)  <  B  <  oo. 

Considering  (170),  the  0(f_1)  term  in  the  asymptotic  expansion  is 

1, 
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(182) 


(183) 


-[S~{K,(z,y)(S+fh(t,2,y)  -  4o(M))}]  =  0J 
Again  we  introduce  the  corrector  Xv(z,y),  and  take  t)i  in  the  form 

»h(t.*.y)  =  X*(z.v)4o(M) 

which  gives  the  equation  for  the  corrector 

5'{A',(«,y)[5+xn(*,y)  -  1]}  =  0 

or 


(184) 


(185) 


(186) 

(187) 


5  {K.{z,y)S+xn(z,y)}  =  Kt{z,y)  -  K,[z,y  -  1) 

By  hypothesis  the  right  side  in  (187)  is  periodic  in  y  and  has  zero  average  (182).  Hence, 
(187)  has  a  periodic  solution,  unique  to  an  additive  constant. 

Continuing,  the  O(c0)  term  in  (169)  is 


S+{rKi(z,y)S  <f>2(t,z,y )}  +  Kt(z,y){S+r}l(t,  z,y)  -  <fii(t,z,y)] 


+  Kt(z,y)~(t,z)  +  S+{rKt(z,y)^-<t>  i(t,z,y)} 


(188) 


+S+{rKt(z,y)^-i4>0(t,z,y)}  +  -^{rKb[z,y  +  l)}-^-<t>0{t,z) 


d 2 


{rKt{z,y  +  l)}^o(M)  ~ 


aVo 


=  0. 


'  Qz2  K  —  •  -UTWV-.-;  Qt2 

This  should  be  regarded  as  an  equation  for  <f>2  as  a  function  of  y  with  (t,z)  as  parameters. 
In  this  sense  the  solvability  condition  is  as  before,  the  average  of  the  sum  of  all  terms  on 
the  left  in  (188),  except  the  first,  should  be  zero.  We  must  choose  <j>o  so  that  this  in  fact 
occurs;  and  that  defines  the  limiting  system. 

Using  the  correctors  (179)  (185),  we  must  have 


Average  (S/K^5  -  ^^[S+(rJC,(r,i/))  +  S*{rK.{z,y)x,(z,l/))} 


(189) 
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-4>Q[^L(rKt(z,y  +  1))  -f  S+(riiCk(z,y)^x#(«,y)) 

+-K'.(«,y)(5+x^(^,y)  -  x*(*»y))]}  =  o 

Defining  the  functions  EI{z),G(z)  by  the  associated  averages  in  (189),  the  averaged  equa¬ 
tion  is 

^£  =  l<£/w^)  +  GW^--ffWo  (190) 

which  is  the  angular  component  of  the  Timoshenko  beam  system  (Crandall  et  al.  [27]  p. 
348). 

Arguing  in  a  similar  fashion,  we  can  derive  the  equation  for  the  macroscopic  approxi¬ 
mation  displacement  of  the  lattice  system  in  terms  of  the  “equivalent”  displacement  r}0(t,z) 
in  the  Timoshenko  beam  system 

^  =  <191> 


4.2.3  Summary 

We  have  shown  that  a  simplified  model  of  the  dynamics  of  the  truss  with  rigid  cross 
sectional  area  may  be  well  approximated  by  the  Timoshenko  beam  model  in  the  limit  as 
the  number  of  cells  (proportional  to  L/l)  becomes  large.  The  continuum  beam  model 
emerges  naturally  in  the  analysis,  as  a  consequence  of  the  periodicity  and  the  scaling. 

To  compute  the  approximate  continuum  model,  one  must  solve  (180)  and  (187)  (numer¬ 
ically)  for  the  correctors  and  then  compute  the  parameters  in  (190)  (191)  by  numerically 
averaging  the  quantities  in  (189)  (and  its  analog  for  (170))  which  involve  the  correctors 
and  the  data  of  the  problem. 

4.3  Homogenization  and  Stabilizing  Control  of  Lattice  Struc¬ 
tures 

In  this  subsection  we  show  that  the  process  of  deriving  effective  “continuum”  approxi¬ 
mations  to  complex  systems  may  be  developed  in  the  context  of  optimal  control  designs 
for  those  systems.  This  procedure  is  more  effective  than  the  procedure  of  first  deriving 
homogeneous  -  continuum  approximations  for  the  structure,  designing  a  control  algorithm 
for  the  idealized  model,  and  then  adapting  the  algorithm  to  the  physical  model.  In  fact, 
separation  of  optimization  and  asymptotic  analysis  can  lead  to  incorrect  algorithms  or 
ineffective  approximations,  particularly  in  control  problems  where  nonlinear  analysis  (e.g., 
of  the  Bellman  dynamic  programming  equation)  is  required. 

We  shall  apply  the  combined  homogenization  -  optimization  procedure  described  in 
subsection  4.1  (based  on  (Bensoussan,  Boccardo,  and  Murat  [59]))  to  the  problem  of  con¬ 
trolling  the  dynamics  of  lattice  structures  like  the  truss  structure  analyzed  in  the  previous 
subsection.  We  shall  only  formulate  a  prototype  problem  of  this  type  and  discuss  its 
essential  features. 
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Figure  3:  Truss  with  transverse  actuator  forces. 


f- 


i 

■j* 


m 


Consider  the  model  for  the  lattice  structure  analyzed  in  subsection  4.3  with  control 
actuators  added.  The  truss  shown  in  Figure  1  is  again  constrained  to  move  in  the  plane 
and  torsional  motion  is  excluded  to  simplify  the  model  and  confine  attention  to  the  basic 
ideas.  Now,  however,  we  include  a  finite  number  of  actuators  acting  to  cause  transverse 
motions.  The  truss  with  actuator  forces  indicated  by  arrows  is  shown  in  Figure  3.  The 
corresponding  discrete  element  model  is  shown  in  Figure  4. 

Suppose  that  the  physical  actuators  act  along  the  local  normal  to  the  truss  midline 
as  shown  in  the  figures,  and  that  the  forces  are  small  so  that  linear  approximations  to 
transcendental  functions  (e.g.,  sin<f>i  —  etc.)  are  valid.  Then  the  controlled  equations 
of  motion  of  the  discrete  element  system  are  (recall  equation  (164)) 


r 


d'4>\ 
dt 2 


-  *;}  +  v+{jr;v<-tf(()> 


(192) 


S  =  -«(<)]) +  £«(*,  •>,(<) 

;=i 

where  the  notation  in  (165)  has  been  used, 

*(«.»=  {J  \t3j  (193) 

and  ij,j  =  l,...,m  are  the  locations  of  the  actuators.  Hence,  if  5(i,t;)  =  0  for  all 
j  —  1, . . . ,  m  there  is  no  actuator  located  at  the  i,h  point  which  corresponds  to  the  physical 
point  x  €  [0,  L).  The  number  m  of  actuators  is  given  at  the  outset  and  does  not,  of  course, 
vary  with  the  scaling.  The  control  problem  is  to  select  the  actuator  forces  as  functions  of 
the  displacements  and  velocities  of  components  of  the  structure  to  damp  out  motions  of 
the  structure.  Measurements  would  typically  be  available  from  a  finite  number  of  sensors 
located  along  the  structure.  We  shall  not  elaborate  on  this  component  of  the  model,  and 
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Figure  4:  Discrete  element  model  of  the  controlled  truss. 

shall  instead  assume  that  the  entire  state  can  be  measured.  To  achieve  the  stabilization, 
we  shall  associate  a  cost  functional  with  the  system  (102).  Let 

«(0  =  K(0.-»u*(<)]T  (194) 

be  the  vector  of  control  forces,  and 

^h-)]=  riii+imr+tiW)? 

j°  i=i 

+  +  Al#(01*  (195) 

+  £*(*»*/)«*  W}*-***1 

j=i 

where  (a,-,  6,)  and  (a,, /9.)  are  non-negative  weights.  Formally,  the  control  problem  is  to 
select  S(i,  i})u}(t),  i  =  1, . . . , N,j  =  1, . . . ,  m  to  achieve 

«(•)]  (1»6) 

subject  to  (192)  (193)  and  the  appropriate  boundary  conditions.  The  case  7  — ►  0  corre¬ 
sponds  to  stabilization  by  feedback. 

The  analysis  of  this  control  problem  is  based  on  the  scaling  used  in  subsection  4.3, 
equations  (164)  -  (171).  Let  r  =  t/e  be  the  fast  time  scale,  then 

.>>(•)]  =  rcEWAwr  +  MWl* 

Jo  i= 1 

+  fl,**(i,(r)]*  +  Ac1|«;(r)]* 


inf 

«()  1 


(107) 
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+  EHi>i>W(r)}*-“'Tdr 

i= * 

with  #(r)  =  $(tr),  etc. 

Let  {<f>,4>,ri,ri)  be  the  state  vector  of  the  system  (102)  with  4>  =  [^i, . . .  ,^jv]T  and 
similarly  for  the  other  terms.  Let  V  =  Vt<"l(4>y4>,r},rj)  be  the  optimal  value  function  for 
the  problem  (102)  (107).  Then  the  Bellman  equation  associated  with  (102)  (107)  is 

<D^  +  W.l 


t=l  r  r 

+  <i;{-v  |*;(  v-.,.  -  *.)|}  v„, 

«=1 

U>€U^  j=i 

+t  +  M?  +  e2 3{on<f>2  +  0,f)2)}  -  nV  =  0. 

•=1 


Remarks: 


(198) 


1.  Note  that  the  minimization  in  (108)  is  well  defined  if  the  admissible  range  of  the 
control  forces  is  convex  since  the  performance  measure  has  been  assumed  to  be 
quadratic  in  the  control  variables  6(t,  tj)uj. 

2.  Since  we  have  not  included  the  effects  of  noise  in  the  model,  the  state  equations  are 
deterministic  and  the  Bellman  equation  (108)  is  a  first  order  system.  To  “regularize” 
the  analysis,  at  least  along  the  lines  followed  in  conventional  homogenization  analysis, 
it  is  useful  to  include  the  effects  of  noise  in  the  model  and  exploit  the  resulting 
coercivity  properties  in  the  asymptotic  analysis. 

3.  If  we  introduce  the  macroscopic  spatial  scale  z  =  ci,  the  mesh  {x<},  and  the  variables 

=  tf(0.  etc.  (100) 

then  the  sums  may  be  regarded  as  Riemann  approximations  to  integrals  over  the 
macroscopic  spatial  scale  z.  The  asymptotic  analysis  of  (108)  with  this  interpreta¬ 
tion  defines  the  mathematical  problem  constituting  simultaneous  homogenization  - 
optimization  for  this  case. 
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4.4  Effective  conductivity  of  a  periodic  lattice 

In  this  subsection  we  consider  a  version  of  a  heat  conduction  problem  treated  by  Kunne- 
mann.  Simple  expressions  for  thermal  properties  of  composite  materials,  have  been  derived 
in  the  past  using  homogenization  techniques.  The  derivation  of  effective  conductivities  for 
discrete  structures  is  useful  for  assessing  the  behavior  of  such  structures  in  variable  envi¬ 
ronmental  conditions. 


4.4.1  Problem  definition 

Let  Z  —  {0,  ±1,±2, ...}  and  Zd  =  Z  x  ...  x  Z  ( d  times)  be  a  d-dimensional  lattice.  Let 
e  >  0  be  a  number  small  relative  to  1.  We  want  to  describe  the  effective  conduction  of 
thermal  energy  on  the  e-spaced  lattice  eZd.  Let  e,  =  (0,0, . .  .0, 1,0, . .  .0)r  with  1  in  the 
ith  position,  »  =  1,2,  . . .  ,d.  If  z  is  a  point  in  eZd,  then  x  ±  ee,,  1  <  »  <  d,  are  the  nearest 
neighbors  of  x.  Let  a±(x),x  6  Zd,  1  <  i  <  d,  be  the  two  functions  defined  on  the  lattice, 
and  assume 


Next  let 


o,-(x)  =f  a+(x)  =  a,_(x  +  e<),x  e  Zd,l<i  <d 
0  <  A  <  c,(x)  <  B  <  oo,  Vx  E  Zdt  1  <  *  <  d 
a,(x)  is  periodic  with  period  i  >  1 
in  each  direction,  1  <  *  <  d. 

®<±(x)  =  <*»±(-)>x  €  tZdy  1  <i  <  d. 


(200) 

(201) 

(202) 

€ 

(203) 


Equation  (201)  means  that  the  conduction  process  is  reversible  and  that  the  conduc¬ 
tivity  d(x)  is  a  “bond  conductivity,”  i.e.,  independent  of  the  direction  in  which  the  bond 
(x,x  +  e,)  is  used  by  the  process.  Equation  (203)  means  that  the  configuration  of  bond 
conductivities  ai±  (•)  on  tZd  is  simply  a ,±  (•)  on  tZd  “viewed  from  a  distance.”  Assumption 
(202)  imposes  a  regularity  condition  on  the  physics  of  the  conduction  process.  An  assump¬ 
tion  like  this  is  essential  for  existence  of  a  limit  as  e  — ►  0.  In  one  dimension  the  situation 
is  illustrated  in  Figure  5  and  Figure  6.  A  system  similar  to  this  with  random  bond  con¬ 
ductivities  was  treated  by  Kunnemann  [78]  by  imposing  some  ergodicity  properties  on  the 
bond  conductivities. 

One  can  associate  with  this  system  a  random  (jump)  process 

{X'{t,x),t  >  0, x  6  eZd} 

on  the  e-spaced  lattice.7  In  effect,  as  e  —*  0  ,  {A-*}  converges  to  a  Brownian  motion  on 
the  lattice;  and  the  main  result  of  the  analysis  is  an  expression  for  the  diffusion  matrix 
Q  =  =  1,2  of  this  process.  This  matrix  describes  the  macroscopic  diffusion 

of  thermal  energy  in  the  system.  It  is  the  effective  conductivity. 

6The  period  may  be  different  in  different  direction*. 

T Definition  of  thia  procese  is  not  necessary  for  the  analysis,  but  it  bolsters  the  intuition. 


.  /Vv  oW/  V  v  v  ‘s.-r. v.  v _y jr.  f. *  -f 
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Figure  5:  Conductivity  on  unsealed  lattice  with  period  1  =  6. 


Figure  6:  Conductivity  on  e-scaled  lattice,  y  =  ex,  x  €  Z,  period  U  =  6e. 
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We  shall  carry  out  the  asymptotic  analysis  of  this  system  in  the  limit  as  c  — ►  0  using 
homogenization.  Let 


(V;  «)(*)  =  -[u(i-c«0  -  u(z)] 

(VJ+.)(*)«i[«(«  +  €^) -•(*)! 


(204) 


xeeZd,  1  <  *  <  d, 

for  any  u  square  summable  on  eZd  or  square  integrable  on  kd  with  e<  the  tth  natural  basis 
vector  in  R.d.  Then 

^  =  rMf)v;M«.*))  ;  poo) 


(206) 


=  Llu‘(t,x) 

is  the  diffusion  equation  on  the  e-spaced  lattice  with  density  u'(t,  i)  and  conductivity 
tii(x/e).  We  are  interested  in  an  effective  parameter  representation  of  the  thermal  conduc¬ 
tion  process  as  c  — *  0. 

Remark:  Although  probabilistic  methods  are  not  required  in  the  analysis,  the  associ¬ 
ated  probabilistic  framework  has  a  great  deal  of  intuitive  appeal.  The  operator  IS  may  be 
identified  as  the  infinitesimal  generator  of  a  pure  jump  process  Xc(«)  in  the  “slow"  time 
scale  a  =f  e2t;  (Breiman  [63]).  Moreover,  Ll  is  selfadjoint  on  tZd  with  the  inner  product 

(/,?)  =f  53  /(*)ff(*)-  (207) 

z€Zd 

Hence,  the  backward  and  forward  equations  for  the  process  X'(s)  are,  respectively, 

— I'-1*-1  =  [iV(»,l|)IW  (208) 

.  IW, <!*))(*) 

So  the  process  is  “symmetric”  in  the  sense  of  Markov  processes  (Breiman  [63]). 

The  asymptotic  analysis  of  (206),  when  interpreted  in  this  context,  means  that  as  the 
bond  lattice  is  contracted  by  e  and  time  is  sped  up  by  t-2,  the  jump  process  {^(s)} 
approaches  a  diffusion  process  with  diffusion  matrix  Q.  In  other  words,  on  the  microscopic 
scale  thermal  energy  is  transmitted  through  the  lattice  by  a  jump  process;  but  when  viewed 
on  a  macroscopic  scale  the  energy  appears  to  diffuse  throughout  the  lattice.  The  micro¬ 
scopic  physics  are  described  in  (Kirkpatrick  [66])  and  (Kittel  [67]).  The  approximation 
developed  below  for  a  periodic  lattice  is  similar  to  the  one  developed  by  Kunnemann  for  a 
random  lattice.  This  similarity  demonstrates  the  robustness  of  the  method,  and  the  limted 
dependence  of  the  macroscopic  properties  of  the  medium  on  the  details  of  the  microscopic 
variations  of  the  structure. 

Because  the  basic  problem  (206)  is  “parabolic,”  we  can  introduce  the  probabilistic 
mechanism  and  make  use  of  it  in  the  analysis.  In  the  “hyperbolic,”  structural  mechanical 
problems  we  treated  before  this  device  is  not  available. 
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4.4.2  Asymptotic  analysis-homogenisation 

The  essential  mathematical  step  is  to  show  strong  convergence  of  the  semigroup  of  L‘,  say 


r(f)  =f  eLU— >T{t)  =  tu  as  e 


4if  mu 


and  to  identify  the  limiting  operator 


L  ' 

This  is  accomplished  by  proving  convergence  of  the  resolvents 


for  a  >  0,  [— U  +  a 


[— L  +  a]  1  as  £  —*•  0 


That  is,  if  /  is  a  given  function  and 


«'(•)  =  | -V  +  a]-1/ 

«(•}  d=  [-1  +  c)-'f 


(209) 


(210) 


(211) 


(212) 


then  u*  — ♦  u  in  an  appropiate  sense. 

The  method  of  multiple  scales  will  be  used  to  compute  the  limit.  Because  the  conduc¬ 
tivities  Oj(x)  in  (206)  do  not  depend  on  time,  we  may  work  directly  with  U  rather  than 
the  parabolic  PDE  (206)  (cf.  (Bensoussan,  Lions,  and  Papanicolaou  [60])  Remark  1.6,  p. 
242).  The  method  of  multiple  scales  is  convenient  because  it  is  a  systematic  way  of  arriving 
at  the  “right  answers”  -  something  which  is  not  always  simple  in  this  analysis. 

Bearing  in  mind  (212),  we  consider 


(LV)(z)  =  /(x) 


with  u*(x)  in  the  form 


j  jr  jj 

1i*(x)  =  Uq(x,  -)  +  £Ui(x,  -)  +  £JU2(X,  -)  +  ••. 


(213) 


(214) 


with  the  functions  «>(x,  y)  periodic  in  y  £  tZ*  for  every  j  —  0, 1, . . ..  (As  it  turns  out  the 
boundary  conditions  are  somewhat  irrelevant  to  the  construction  of  “right  answers.”)  To 
present  the  computations  in  a  simple  form,  it  is  convenient  to  introduce  y  =  x/£,  to  treat 
x  and  y  as  independent  variables,  and  to  replace  y  by  x/f  at  the  end. 

Recall  the  operators  VJ*  from  (206).  Applied  to  a  smooth  function  u  =  u(x,  x/e),  we 
have 

(VJ"u)(x,y)  =  ^[u(x-  £e„y-  e.)  -  u(x,y)|  (215) 

=  ^(u(x,y  -  e.)  -  u(x,y)]  +  j[u(x  -  £e„y  -  e.)  -  u(x,y  -  e,)] 

1,__  w  v  du ,  .  I5su.  , 

=  -y,V-  e<)  +  e2dx*(X,V  ~  ^  +  ) 


•>  VV  >  W  V  V  *,•  V  V  V 
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8 

s 

t 

5 

*  ta 


where  on  functions  <t>  —  4{y) 


Defining 
we  also  have 


(V,  ^)(v)  =  ^(y  - «.)  -  <t>(v) 
(vi»(y)  =  ttv  +  «i)  -  4>{y) 

(VJ+u)(x,y)  =  ^(V,tu)(z,y)  +  ^-(x.y +  e.) 

+^!^(*.V  +  «.)  +  0(e2}.  ; 


(216) 

(217) 

(218) 


Now  we  substitute  (214)  into  (213)  and  use  the  rules  (216)  (217).  Equating  coefficients 
of  like  powers  of  e,  this  leads  to  a  sequence  of  equations  for  uq,ui,  . . ..  Specifically,  (using 
the  summation  convention) 


fa 


i 

§ 

* 

V 


(LV)(x,y)  =  -Vr[«.(y)V;+u‘j 
=  ^v,:[fl<(y)v,t«o(i,y)]  -  ^'K(y)^:(x,y  +  «)] 
-  ^cVJ-[o<(y)^(x,y  +  «<))  +  °(0 
-jVf[a,(y)Vtu1(x,y)] 
-eVJ-[a,(y)|^(x,y  +  e,)]  +  0(0 

-vr(«.(y)v.ttt»(*.y)]  +  o(0  =  /(*) 

That  is,  labeling  each  term  by  its  order  in  e 


(219) 


(«"*)  Vr[a,(y)Vtu0]  =  0 


(220) 


(e'1)  cVJ  la,(y)|^(x,y  +  e,)J  +  [a,(y)Vtu,(x,y)J  =  0  (221) 

and  (recall  eVJ*  is  0(1)  in  0 


& 


(c°)  icVJ  K(y)^(x,y  +  e,)]-£V,‘  [a,(y)^-(x,y  +  e<)) 


(222) 


>2 

£ 


* 


£ 

v 


-  v.r  («» ( v)  <*2  » y )  ]  =  /(*) 

From  (220)  we  have 

«<(y  -  «<)ltto(*.y)  -  t*o(x,y  -  «0) 
-Oi(y)[u0(z,y  +  e,)  -  u0(x,y)j  =  0 


(223) 


3? 
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If  we  take  u0(x,y)  =  u0(i),  this  is  trivially  true;  and  (221)  simplifies  to 

«vrk(v)|^(x)l  +  v<"  [«<(»)  v<+«i(*ty)J  =  °-  (224) 

At  this  point  we  introduce  “correctors.”  That  is,  we  assume 

«x  (*,  y)  =  Y,  x*  (y)  jrr  +  (z)  (225) 

*=i 

with  x*(*)  the  correctors.  Using  this  in  (224),  we  have  (again  using  the  summation  con¬ 
vention) 

vrMvW x*(y)]|^  +  My  -  **)  -  °*(y)]|^  =  0  (226) 

If  we  take  x*(y)  ®s  the  solution  of 

vrMyW x*(y)|  +  (a*(y  -  tk)  -  My)]  =  o  (227) 


(we  have  to  verify  the  well-posedness  of  (227)),  then  (226)  is  satisfied.  (The  term  ui(x)  is 
determined  (formally)  from  the  O(e)  term  in  the  system  (214)  (219).) 

Regarding  the  well-posedness  of  (227),  note  that 


vrMy)Vt*(y)]  =  v>(y)  (228) 

has  a  periodic  solution  on  eZ  which  is  unique  up  to  an  additive  constant  iff  the  average  of 
the  function  t/>(y)  over  a  period  (e£)  is  zero;  i.e., 


=  yX]^(y  +  M  =°>n  =  l,2,...,d. 

L  k=i 


(229) 


This  condition  clearly  holds  in  (227),  and  so,  X*(y)  is  well  defined  (up  to  an  additive 
constant) . 

We  shall  determine  the  equation  for  u0(x)  by  using  (225)  (227)  in  (222).  Using  the 
Kronecker  delta  function  we  have 


icvrwsOM 


d2u0 

dxidxk 


eeV;  (ai(y)xk(y +  «<)] 


3Ju0 

dxidxk 


f(x) 


=  -  Vr[o.(»)V*Xl]}g^l-  (230) 

-  V,'[a,(v)V*u,l  =  f(x) 

The  term  in  braces  is  zero  from  (227).  To  obtain  the  solvability  condition  (229)  for  uj  in 
(230),  we  introduce  the  average 


5  9.*  =  symmetric  part  {  —  V ~  (a,  (y) Xk (y) ]  > 


(231) 
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Then  solvability  of  (230)  for  U2  gives  the  equation 


1  ^  d3u0  ,,  s 


(232) 


And  this  is  the  diffusion  equation  which  defines  the  limiting  behavior  of  the  system  (213) 
in  the  macroscopic  x-scale  in  the  limit  as  c  — ►  0. 


We  can  justify  the  asymptotic  analysis  by  using  energy  estimates  or  probabilistic  meth¬ 
ods  as  in  (Bensoussan,  Lions,  and  Papanicolaou  [60]).  (See  also  Kunnemann  lobal  set  key 
[78]).)  We  shall  omit  this  analysis  here. 


4.4.3  Summary 


Returning  to  the  original  problem  (206)  for  the  evolution  of  thermal  energy  on  a  micro¬ 
scopic  scale,  we  have  shown  that  the  thermal  density  u‘(f,x)  — »  u0(f,x)  as  f  — ►  0  (in  an 
appropriate  norm)  where 

(2*3) 


2  1  dxidxj 


9*i  =  “T  !£WMy)xh(v)]  +  Vk  [fl,(y)x.(y)]> 

Lk= i 


(234) 


with  the  correctors  Xk» k  =  1,2, ...  ,d,  given  by 


J2  MyW x<(y)l  =  ~Mv  -  **)  -  a*(y)] 


(235) 


k  =  1,2  ,...,d 


To  compute  the  limiting  “homogenized”  model  (233),  one  must  solve  the  system  (235) 
(numerically)  and  then  evaluate  the  average  (234). 

The  fact  that  the  original  problem  (206)  is  “parabolic”  (i.e.,  it  describes  a  jump  random 
process),  enables  us  to  exploit  the  associated  probabilistic  structure  to  anticipate  and 
structure  the  analysis.  In  this  way  we  can  anticipate  that  the  limit  problem  will  involve 
a  diffusion  process.  In  fact,  the  arguments  used  are  entirely  analytical  and  the  limiting 
diffusion  (233)  is  constructed  in  a  systematic  way.  It  is  not  postulated.* 


'Probabilistic  arguments  can  be  used  (Bensoussan,  Lions,  and  Papanicolaou  (GO],  Chapter  S);  and  they 
have  some  advantages. 
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Figure  1:  SCOLE  Model 


5  Examples 

5.1  Continuum  modeling  for  SCOLE  problem 

The  SCOLE  design  challenge  provides  a  useful  baseline  for  comparison  of  continuum  mod¬ 
eling  and  control  design  for  vibration  problems  typically  encountered  in  large  spacecraft. 
The  configuration  is  discussed  in  detail  in  [79]  and  consists  of  a  large  reflector  antenna 
mounted  on  the  end  of  a  long  truss  which  is  attached  to  standard  pallet  in  the  shuttle  or- 
biter  cargo  bay.  The  simplified  model  discussed  in  [84]  includes  a  standard  Bernoulli-Euler 
beam  model  for  the  flexible  truss  and  both  the  space  shuttle  and  the  antenna  are  modeled 
as  rigid  bodies.  Thus  the  problem  can  be  considered  as  a  “barbell”  with  flexible  bar  con¬ 
necting  the  masses.  In  this  section  we  provide  a  linear  hybrid  model  for  this  configuration 
by  considering  motion  in  a  plane.  The  configuration  is  shown  in  Figure  5.1. 

The  required  equations  of  motion  can  be  easily  derived  from  by  application  of  Hamil¬ 
ton’s  principle.  In  terms  of  the  dynamic  variables  listed  in  Table  2  we  can  write  the 
Lagranian  as 


•1 


* 


SEI-TR-86-14 

torque  applied  to  shuttle 
angular  attitude  of  shuttle 
effective  shear  modulus  of  truss 
mass  density  of  truss 
cross  sectional  area  of  truss 
modulus  of  elasticity  of  truss 
moment  of  inertia  of  truss 
truss  length 

mass  and  moment  of  inertia  for  antenna 
mass  and  moment  of  inertia  of  shuttle  (about  composite  CG) 
distance  from  point  of  truss  attachment  to  CG 
longitudinal  coordinate  along  truss 
lateral  displacement  of  truss 
angular  rotation  of  cross  section  of  truss 

Table  2:  Notation  for  SC  OLE  planar  motion  model 

+  ijWi.  (^  +  (i  +  +  • 

Then  by  application  of  the  standard  variational  calculus  [27]  to  Lagrange's  equation, 

i_  (dL\  _ 
dt  \dq  J  dq  ’ 


the  following  equations  of  motion  result.  The  Distributed  Parameter  Subsystem  (DPS)  is 

d2r]  kG  ( d2r) 

dt 2  p  \dz2 

d2<i>  E  d2<t> 

dt2  ~  p  dz2  + 

-0 -<*+*>» 

(237) 

with  boundary  conditions 

*?M)  =0, 

0)  =  0 

(238) 

r)(t,L)  =  r?i(t), 

4>{t,L)  =  <j>L. 

The  Lumped  Parameter  Subsystem  (LPS)  is 

Id  -f  Mifji  +  Ji4>i  = 

•“(SI..-SIJ 

(230) 

(240) 

T 

6 

kG 

P 

A 
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MJ 
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V(z) 

4>(z) 
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Figure  2:  Subsystem  interaction  for  SCOLE  linear  model 


Ml(L  +  R)6  +  Mltil  ~  kGA4>l  =  -kGA  — ^ 

oz 

-  »  dd> 

hO  +  h<t>  =  -  ~ 


(241) 

(242) 


The  Distributed  Parameter  System  (DPS)  interacts  with  the  Lumped  Parameter  Sys¬ 
tem  (LPS)  at  the  boundaries  as  shown  in  Figure  5.1.  For  the  purposes  of  illustration  we 
consider  simplifying  this  model  by  assuming  the  shuttle  mass  is  very  much  larger  than  that 
of  the  antenna  reflector. 


5.2  Cantilevered  Beam  Models 

A  natural  simplification  of  the  above  model  is  to  assume  the  shuttle  body  is  so  massive  with 
respect  to  the  antenna  reflector  that  the  vibration  dynamics  can  be  effectively  studied  by 
assuming  the  shuttle  is  rigidly  fixed  in  space.  In  this  section  we  consider  the  ‘cantilevered’ 
vibration  model  for  a  flexible  beam  with  a  tip  mass.  This  is  a  standard  problem  which  is 
derived  in  detail  in  [27,  pp.  348].  For  control  of  cantilevered  vibration  modes  we  include  a 
control  force  at  the  antenna  relector  fu(t)  which  acts  orthogonal  to  the  longitudinal  axis 
of  the  truss  in  it  quiescent  state.  We  state  the  equations  of  motion  as: 

Distributed  Subsystem 


.d2f] 
pA  dO 


-  £b(§H] 


d  „  Td<f> 
—  EI-~ 
dz  dz 


+  kGA 


ft-*) 


(243) 

(244) 
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Lumped  Subsystem  (at  z  =  L) 
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Mx.r?L  +  xGA  -  «^>l  j  =  fu 


(0 


El^jr  +  hh  =  0 

dz 


(245) 

(246) 


and  with  boundary  conditions  at  z  =  0 

t?(4,0)=0  #(4,0)  =0. 


(247) 


As  in  Section  3  we  motivate  the  Euler-Bernoulli  model  by  introducing  the  limiting 
argument  leading  to  the  equations 


d2 


_ 

pA  dt 2  5z2 


(*'B) 


with  boundary  conditions  at  x  =  L 


d^n 

Mlth  =  EI-q^  +  /u(0> 


£/S?  =  0 


dz1 


at  x  =  0 


r](t,0)  =  0 


di  -  ° 


To  illustrate  the  modeling  techniques  for  hybrid  systems  and  the  computation  of  ef¬ 
fective  state  space  models  we  consider  the  ‘long,  thin  beam’  approximation  discussed  in 
Section  2.  For  this  simple  example  we  will  use  the  approximation  (discussed  in  Section  2 
as  the  ‘string’  model)  in  which  the  Timeshenko  equations  are  reduced  by  neglecting  the 
rotation  angle  of  cross-section  <f>  =  0  in  (54).  The  model  for  the  hybrid  system  consists  of 
equation  (69)  for  0  <  z  <  L,  subject  to  boundary  conditions  at  z  =  L, 


mL 


*^  +  KGA?>ML  =  m 


dt 2 


(248) 


and  at  z  =  0, 


,M)  =  =  o. 


dz 


(249) 


Now  (69)  can  be  written  in  form  (100)-(99)  by  a  particular  choice  of  distributed  state: 
Xi(t,z)  =  |7,»7]t  where  (69)  becomes  with  o2  =  kG/p 


dq(M) 

dt 

di{t,z) 


=  a 


drr(M) 


dt 


=  Q 


dz 

dr]{t,z) 


dz 


(250) 
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7r(s,.Z,U>)  =  |  ^ 


RIGHT > 


0  <  Z  <  W 

W  <  Z  <  L 


Hbc{$,z)  = 


cosh  ^ 
sinh  — 

_ Cr 

sinh  — 

a 


The  final  hybrid  model  can  be  written  in  terms  of  (106)  by  identifying  the  following 
terms;  from  (109),  the  hybrid  Green’s  function  Gr  can  be  expressed  using  (109)  given  the 
terms 


-HtQ,P  = 


—s  sinh  ^  s  cosh  ^ 
-  sinh  bb.  cosh  bb. 

_ a _ a 


cosh  —  —  mis2  sinh  — 

a  u  a 

—  sinh  ~  (w  —  z)  —  sinh  £  (tu  +  z)  cosh  £  (w  —  z)  +  cosh  £  (tu  +  z) 

-  _  cosh  £(tu  —  z)  —  cosh  i(u>  +  2)  —  sinh  £(tu  —  z)  +  sinh  £(u>  +  z) 

2  2  sinh  ^  (cosh  ^  —  m/,s2  sinh  ^) 

From  (108),  we  can  write 


R{s,z)  = 


—mrs  sinh  — 

r° 

-m>  sinh  — 

°L 

-mi  cosh  ~ 
—mL  sinh  ^ 
cosh  —  —  : 


—  cosh  — 

a 

—mLs  cosh^ 
—mis  cosh  — 
~mLs  cosh^ 


mis1  sinh  ^ 


Finally,  from  (111) 


—5  cosh  — 
_  it 


—  cosh  — 

Cl 

—  mis2  cosh  — 

~  —mis2  sinh  ^ 

’  cosh  —  -  mx,sJ  sinh  — 

The  above  calculations  were  carried  out  using  SMP.  Some  considerable  algebraic  re¬ 
duction  was  necessary  to  achieve  the  final  forms. 


5.3  Control  Design  for  Simply  Supported  Beam 

For  illustrative  purpose  the  control  design  for  a  simply  supportted  (‘pinned-pinned’)  beam 
is  considered  in  some  detail.  Continuum  modeling  begins  with  standard  Bernoulli-Euler 
beam  equation.  We  develop  a  canonical,  first-order  form  and  identify  suitable  state  space 
coordinates  and  identify  boundary  conditions  which  places  the  model  in  the  canonical  form 
discussed  in  Section  3.  Computation  of  the  required  frequency  response  requires  numer¬ 
ical  evaluation  of  various  transcendental  terms.  Numerical  aspects  of  this  are  discussed. 


» 


.vlvS.’/.'.-I' 


SEI-TR-86-14 


61 


Finally,  a  distributed  control  is  computed  by  the  method  described  in  Section  2.  We  high¬ 
light  the  computational  aspects  of  the  procedure  and,  in  particular,  the  use  of  a  symbolic 
manipulation,  computer  algebra  system  to  support  the  modeling  and  control  design. 

Testing  of  the  distributed  control  law  proceeds  in  the  frequency  domain.  This  is  an 
advantage  of  the  approach  since  frequency  response  data  is  directly  available  as  part  of 
the  modeling  and  control  computations.  The  ultimate  concern  for  design  of  active  control 
of  elastic  structures  in  space  is  to  provide  adequate  stability  margins  in  the  face  of  the 
inevitable  model  uncertainty.  The  most  standard  methods  for  testing  stability  and  assess¬ 
ing  the  robustness  of  control  laws  involve  frequency  domain  tests.  Such  tests  (based  on 
Nyquist  stability  theory)  are  firmly  based  in  engineering  practice  and  proven  reliable.  The 
inherent  robustness  properties  of  the  standard  formulation  for  linear,  quadratic  optimal 
control  which  rest  on  frequency  domain  properties  of  the  control  law  can  be  readily  seen  in 
this  framework.  It  is  essential  for  control  of  distributed  parameter  systems  where  computa¬ 
tions  will  inevitably  involve  approximation  that  any  degradation  of  the  resulting  stability 
margins  be  ascertainable  in  advance  for  the  control  configuration  including  sensors  and 
actuators. 

5.3.1  Continuum  Modeling  and  State  Coordinate  Selection 

The  ‘pinned-pinned’  beam  is  shown  in  Figure  5.3.1.  The  beam  is  modeled  by  a  simplified, 
dimensionless  form  of  the  Bernoulli-Euler  equation; 


w  +  £!»  =  o 

dt2  s  dtdz2  dz*  ' 


(254) 


where  the  Chen-Russell  ‘square-root’  damping  is  included  for  material  dissipation  and 
(254)  is  subject  to  boundary  conditions  at  z  =  0 


V{t,  0)  =  0 


(no  lateral  displacement) 


(no  restraining  moment)  and  at  z  =  L  a  control  torque  r  is  applied, 

V  («,L)  =  0 

d7ri  I 

a,’  .  =  r(,) 


(255) 


(256) 


(257) 

(258) 


Writing  (254)  in  the  form 

d2  I  drj  d2r)  \ 

dO  dz2  \2<  dt  dz2  I 

identifies  a  new  state  variable  -j  and  we  write  (254)  in  the  form 


d2l 
( iz2 

d»l  d2v 
dt  dz2 


(250) 
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Figure  4:  Simply  Supported  ‘pinned-pinned’  Beam 
Now  solving  for  17  in  (260)  and  subtracting  from  (259)  gives 


d*i  _  d 2/i  d2rj 
dt  =  “  dz»‘ 


(261) 


With  the  choice  of  state  coordinates  as  x  =  we  can  write  (259)  and  (261)  in  the 

canonical,  first-order  form  as 


d2  (  0  1  \ 

1  dz ■*  \  -1  2f  J  Z‘ 


(262) 


The  boundary  conditions  (256),  (258)  involve  moments  and  can  be  written  in  terms  of  the 
state  coordinates  from  (260)  as 


T((«,0)-2fi7(t,0)  =  0 

Tf(t,I)  -  2f»j(*,I)  =  f  r{a)da; 

Jo 


(263) 

(264) 


i.e.,  as  discussed  in  Section  3  a  generalized  forcing  function  (torque)  is  evident  and  the 
integration  in  (264)  is  essential. 

Finally,  the  equations  (262)-(264)  can  be  written  in  the  canonical  form  (92),  (93) 
identifying  the  matrix  parameters, 


E.  - 


F  =  0, 

H  = 

0,  £ 

[  1 

0  1 

0 

0 

2?  -1 

0  0 

,r,  = 

0 

1 

0 

0 

[  0 

0 

N 

-1 

E,  =  r,  =  0,£>  =  10,0,0,  J  dl}T. 


Numerical  evaluation  of  the  required  frequency  response  data  for  the  model  in  the  form 
(81)  proceeds  from  (95)-(99).  However,  direct  evaluation  can  lead  to  numerical  instability 
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depending  on  the  required  bandwidth.  To  see  how  the  numercial  instability  arises  we  focus 
on  the  transfer  function  from  the  boundary  control  (torque  at  right  hand  end  of  beam)  to 
the  distributed  state; 

x{s,z)  = 

which  consists  of  two  terms;  viz., 

v{s,z)  =  hnT(s,z)f{s), 

7(5,2)  =  h7r(s,z)f(s). 

An  independent  analysis  shows  that  these  transfer  functions  can  be  written  as 
sinh  \iz  sin  \\L  —  sin  Xyz  sinh  A2L 


^f|r(s,2)  — 


(Af  +  Aj)  sin  AiL  sinh  \2L 


(265) 


.  .  (f  —  *Vl  —  ?2)  sinh  A22sin  AiL  —  (f  +  *V I  —  f2)  sin  Aiz  sinh  XtL  .  . 

=  - (A5  +  Ai).inA1L.iahA1L -  (266) 


with 


=  (-?  +  i\/l  -  ?2)s, 

=  (?  +  »Vl  -  ?2)«- 


(267) 


Clearly  evaluation  of  these  transfer  functions  at  «  =  0  by  direct  substitution  leads  to 
indeterminant  expression  § . 

At  this  point  a  computer  algebra  system  (SMP  was  used)  is  most  useful  for  analysis. 
The  system  can  answer  the  question  “is  there  a  singularity  in  hnT  at  s  =  0?”  For  s  —*  0 
we  can  also  let  f  — »  0  which  implies  Aj  =  A2.  Then  SMP  can  readily  perform  the  required 
power  series  expansion  about  s  =  0  as  shown.  The  result  clearly  demonstrates  that  the 
singularity  is  removable.  A  similar  analysis  for  h1T  shows  that  this  term  has  a  first-order 
pole  at  s  =  0  which  confirms  that  the  generalized  torque  input  to  this  model  is  required 
for  the  ‘new  state’  coordinate  and  is  not  physically  evident  in  the  lateral  displacement 
coordinate. 

Clearly,  numerical  evaluation  of  the  frequency  response  data  from  this  model  near 
s  =  0  will  be  numerically  sensitive.  We  have  found  SMP  a  (somewhat)  convenient  tool  for 
generation  of  Fortran  code  for  numerical  evaluation  of  these  expressions.  For  evaluation 
near  5  =  0  we  use  SMP  to  generate  the  required  power  series  expansions  and  automatically 
generate  Fortran  expression  for  their  evaluation  in  ‘low  frequency’  regions. 


5.3.2  Evaluation  of  Frequency  Response  Data 

In  this  section  we  demonstrate  a  general  approach  to  evaluation  of  frequency  response  data 
for  such  beam  models  which 

1.  separates  the  integration  associated  with  the  generalized  torque  input  to  the  model 

2.  identifies  a  reduced  set  of  transcendental  terms  which  can  be  analyzed  independently 
for  numerical  stability  and  from  which  all  frequency  response  data;  viz.,  Hbci  Gt, 
can  be  computed. 
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Using  the  symmetry  of  the  boundary  conditions  for  this  problem  a  further  simplification 
results. 

Equation  (95)  can  be  written  for  the  ‘pinned-pinned’  beam  as 


A  = 


0 

-sG-1 


(268) 


where 


-sG'1  = 


2fs  —s 
s  0 


As  discussed  in  Appendix  A  we  find  it  convenient  for  algebraic  computation  of  the  matrix 
exponential  $(s,  z)  =  e*K  to  use  a  Cayley-Hamilton  expansion; 


$(s,z)  =  il>0(s,z)I4  +  rpi{s,z)  A  +  02(s,z)A2  +  0s(*,z)As.  (269) 


The  Faddeeva  algorithm  (see  Appendix  A)  is  used  to  evaluate  the  scalar  terms  yielding, 


V’o 

1>3 


(A£  —  2fs)  cosh  A 2z  +  (Aj  —  2fs)  cos  Ax z 
(A?  +  A|) 

j^(AJ  +  2fs)  sin  Ais  +  j^(A^  -  2 fs)  sinh  A2z 

(A?  +  k) 

cosh  A22  —  COS  AjZ 

(aFTaI) 

£  sinh  A2z  -  j-  sin  A \z 
(Af+A f)  ’ 


(270) 


where  A2 ,  A2  are  given  in  (267 ) .  These  terms  include  all  transcendental  terms  involved  in  the 
computation  of  Hbc,Gt.  Numerical  stability  is  an  issue  which  can  therefore  be  determined 
independently  for  these  expressions.  We  also  feel  that  this  approach  may  suggest  alternate 
methods  for  approximation  of  transfer  functions  by  rational  approximation.  Such  methods 
may  yield  efficient  and  numerically  stable  schemes  for  frequency  response  modeling. 

The  model  equations  (95)-(99)  for  this  case  can  be  written  in  a  simplified  form  by 
exploiting  the  symmetry  of  the  boundary  conditions.  From  (268)  we  see  that 


so  that  using  (269)  we  can  write  (95)  in  partitioned  form 


*(«,*) 


$i  #2 

♦3 


(271) 
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where 


$i(a,z)  =  xp0I2  -  stl>2G~l 

^o{s,z)  +  2$sip2{s,z) 
stp2{s,z) 

and 

$2(s,z)  =  xl>il2  -  sfoG'1 

V>i(5,z)  +  2  $s\l>3(s,z) 
sip3{s,z) 

At  this  point  symbolic  algebra  system  (SMP)  was  used  to  evaluate  the  limits  of  the 
elements  in  these  matrices  as  a  — ►  0.  The  SMP  session  is  reproduced  as  follows: 

SMP  1.6.0 

10-DEC- 1986  10:56:02.05 

\*  the  matrix  phil  is  parametrized  independent  of  s  */ 

#1  [1] : :  philts.z] 

2  1/2  2  1/2 

cl  Cos [cl  s  z]  ♦  c2  Cosh[c2  s  z] 

•0[1]:  « . . 

2  2 
cl  ♦  c2 

1/2  1/2 
Cos [cl  s  z]  -  Cosh[c2  s  z] 
. >. 


-sip2{s,  z) 
xpo{s>z) 


-sil>3{s,z) 
i/>! (s,z)  ■ 


(272) 


(273) 
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1/2  1/2 
-(Cob [cl  b  z]  -  Cosh[c2  b  z] ) 


2  2 

cl  ♦  c2 

1/2  2 
Cos [cl  b  z]  (2zeta  ♦  cl  ) 

1/2  2 
-  Cosh[c2  s  z]  (2zeta  -  c2  ) 


-» 


2  2 
cl  ♦  c2 

#1 [2] : :  tmp:%; 

\*  rs~2  is  substituted  for  s  prior  to  evaluation  of  the  limit  *\ 

#1  [33 : :  tmp:S[tmp,s->rs“2] 

2  2 

cl  Cos [cl  rs  z]  ♦  c2  Cosh[c2  rs  z]  Cos [cl  rs  z]  -  Cosh[c2  rs  z] 


#0[33 


2  2 
cl  ♦  c2 


2  2 
cl  +  c2 


-(Cob [cl  rs  z]  -  Cosh[c2  rs  z]) 


2  2 
cl  ♦  c2 


Cos  [cl  rs  z]  (2zeta  ♦  cl  ) 


-  Cosh[c2  rs  z]  (2zeta  -  c2  ) 


•» 


2  2 
cl  ♦  c2 


#I[4] : :  tmp[l ,1] 


cl  Cos[cl  rs  z]  ♦  c2  Cosh[c2  rs  z] 


•0(4] : 


2  2 
cl  ♦  c2 


a 


r0 

r, 

A 


A 


£ 


.  v  •  /  /  v  V  v  *.\v  •■.•’."V  -v . 

-'■  A-V-S  -V  V«_ 
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\*  find  th«  limit  rs->0  for  the  1.1  element  of  psil  *\ 
#I[5]::  LimtX.rs.O] 


#0[5] :  1 

\*  find  limits  for  each  element  of  psil  *\ 
•I [6]::  tmp[l,2] 


#0[6]  : 


Cos [cl  rs  z]  -  Cosh[c2  rs  z] 


2  2 
cl  ♦  c2 


•I [7] : :  LimtX.rs.O] 
#0 [7]  :  0 


#I[8]:: 

#0 [8]  : 


tmp[2,l] 

-(Cos [cl  rs  z]  -  Cosh[c2  rs  z]) 

2  2 
cl  ♦  c2 


#1 [0] : :  LimtX.rs.O] 


#0[9]  : 
#I[10] : : 

#0[10] : 


0 


tmp[2,2] 


2 


2 


Cos [cl  rs  z]  (2zeta  ♦  cl  )  -  Cosh[c2  rs  z]  (2zeta  -  c2  ) 


2  2 
cl  ♦  c2 


#1 [11] : :  LimtX.rs.O] 
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\*  perform  similar  limiting  analysis  for  psi2  *\ 
#I[12] : :  tmp:S[phi2[s,z] ,s->rs“2] 


2zeta  (cl  Sinh[c2  rs  z]  -  c2  Sinh[cl  rs  z] ) 


#0 [12] :  {{- 


-  cl  (2zeta  -  c2  )  Sinh[c2  rs  z] 


+  c2  (2zeta  ♦  cl  )  Sintcl  rs  z] 


2  2 
cl  c2  rs  (cl  ♦  c2  ) 


Sinh[cl  rs  z]  Sinh[c2  rs  z] 


2  2 
rs  (cl  ♦  c2  ) 

Sinh[cl  rs  z]  Sinh[c2  rs  z] 

( . ) 


2  2 
rs  (cl  +  c2  ) 


-(cl  (2zeta  -  c2  )  Sinh[c2  rs  z] 


-  c2  (2zeta  ♦  cl  )  Sin[cl  rs  z]) 


2  2 
cl  c2  rs  (cl  ♦  c2  ) 
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#1 [13]  : : 


#0 [13]  : 


tmp[l ,  1] 

2zeta  (cl  Slnh[c2  rs  z]  -  c2  Slnh[cl  rs  z]) 

2 

-  cl  (2zeta  -  c2  )  Sinh[c2  rs  z] 
2 

♦  c2  (2zeta  ♦  cl  )  Sin[cl  rs  z] 


2  2 
cl  c2  rs  (cl  ♦  c2  ) 


# I [14]  : : 


#0 [14] : 


Lim[X,rs  ,0] 

2 

cl  c2  z  (2zeta  +  cl  )  -  cl  c2  z  (2zeta  -  c2 


2  2 
cl  c2  (cl  ♦  c2  ) 

\*  factor  the  resulting  limit  to  simplify  *\ 
#1  [16]  : :  Facm 


#0 [16]  :  z 


#1 [16]  : :  tmp  [1 ,2] 


Sinh[cl  rs  z]  Sinh[c2  rs  z] 
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f 

! 


# I C 18] : :  tmp[2 , 1] 


Sinh[cl  rs  z]  Sinh[c2  rs  z] 

( . ) 

cl  c2 


#0 [18] : 


2  2 
rs  (cl  ♦  c2  ) 


#1 [19] : :  Lim[X.r8 ,0] 


#0 [19] :  0 


#1 [20] : :  tmp [2 . 2] 


2  2 
-(cl  (2zeta  -  c2  )  Sinh[c2  re  z]  -  c2  (2zeta  ♦  cl  )  Sin[cl  rs  z]) 


#0 [20] 


2  2 
cl  c2  rs  (cl  ♦  c2  ) 


#I[21] : :  Lim[X.rs .0] 

2  2 
-(-cl  c2  z  (2zeta  ♦  cl  )  ♦  cl  c2  z  (2zeta  -  c2  )) 

#0 [21] :  . 

2  2 
cl  c2  (cl  *  c2  ) 

#1  [22] : :  Fact*] 


#0 [22]  :  z 


#1 [23] : :  <end> 


Now  from  the  symmetry  of  the  boundary  conditions  we  can  write 


E  +  r$  = 


A  0 


A$i  A$2 


(274) 


where 


-(i  *.)• 


It  is  easy  to  see  that 


+  T<M-1  =  ^  ® 

“*J,*iA  *J»A  • 


I 
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From  (96)  we  write 
M(s,  z)  = 


$2(3,2)!  | 
(A/l(3,2),Af2(3,2)] 


A 


0 


where 


Mi(s,z)  =  ($i(s,2)  -  $2(5,2)$21(5>-£')$l(5»-£'))  A 

M2(s,z)  =  $2(s,z)$t1(s,L)A. 


Finally,  from  (275)  write  (99)  as 


HBc{s,z)  =  $2(s,2)$21 3 *(5.-£')A 


0 

1/s 


Substitution  of  (271),  (275)  into  (99)  yields 


Gt  left(s;z,xv)  =  -$2(s,2)$2  1(s,L)$2(s,L  -  ttf) 

GrRIGHT{s;z,w)  =  [$j  (s,  z)  -  $2(s,  z)^ *  («,  -^)*1  («»  L)]  $2^,  ~ ^)  • 


(275) 


(276) 

(277) 

(278) 


Numerical  evaluation  then  proceeds  by  evaluating  the  transcendental  terms  (270). 
Substitution  into  (272)  and  (273)  is  followed  by  2  x  2  matrix  calculations  (276)-(278). 
Figure  5.3.2  displays  3-dimensional  frequency  response  data  for  the  terms  in  HBc ■  The 
integration  in  h1T  is  evident  in  the  response  at  w  =  0. 


5.3.3  Distributed  Control  Computation 

The  optimal  control  problem  described  in  Section  1  is  now  completely  specified  by  the 
choice  of  a  quadratic  objective  which  can  be  alternately  specified  by  the  choice  of  a  linear 
operator  C  :  l/(fl)  — ►  as  in  (7).  Candidates  for  C  are: 

1.  projection  to  a  modal  or  other  ROM,  finite-dimensional  subspace  of  H (fl)  as  discussed 
in  Section  2 

2.  “point”  control  objective  given  by  the  operation 

y{t)  =  /  x(t,z)6(z  -  zo)dz 
J  n 

3.  weighted,  state  cost 

y(0  =  I  x{t,z)w{z)dz 
J  n 

For  this  problem  we  chose  a  point  control  objective  with  the  observation  or  control  point 
at  2  =  L/2.  It  should  be  noted  that  this  choice  does  not  suppress  the  fundamental  difficulty 
of  computations  for  distributed  parameter  models  by  employing  ROM  approximations.  As 
we  shall  see  the  numerical  problems  are  quite  complex  even  after  we  have  assured  ourselves 
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Figure  6:  Nyquist  plot  for  1/F(tu/)  for  beam  control  computation. 

by  various  tests  that  the  numerical  evaluation  of  the  irrational  transfer  functions  involved 
can  be  performed  with  the  required  precision.  With  this  choice  the  transfer  function  is 

G(s)  =  HBc{s,z  =  L/ 2) 

and  we  must  factor 

I  4-  G*(»u;)G(iu;)  =  F'(iu;)F(iu) 

to  find  its  unique,  causal  spectral  factor  F(iu).  Using  the  frequency  response  model 
described  above  and  the  spectral  factorization  algorithm  discussed  in  Section  2  we  compute 
the  quantity  [F*(tu;)]_1  which  is  in  this  case  a  complex  scalar  valued  function  of  w.  The 
result  is  shown  in  Figures  5.3. 3-5. 3. 3  with  500  samples  computed  over  the  frequency  range 
from  1.0  x  10“10-1.0  x  103. 

Two  observations  can  be  made  from  the  resulting  spectral  factor  which  help  to  confirm 
its  validity.  First,  from  basic  principles  of  linear,  quadratic  control  it  is  known  [23,  pp.  71] 
that  the  unique  causal  spectral  factor  is 

F(s)  =  /  +  KoptR(s;  A)B 


and  enjoys  the  property 

||F(tw)|j  >  1  VweSR. 

Clearly  from  Figure  5.3.3  we  see  that  |F-1(— tu/)  |  <  1  consistent  with  this  property.  Second, 
by  choice  of  control  observation  at  z  =  Lj 2  we  see  that  all  symmetric  modes  in  Hgc  do 
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k-mode 

frequency  (rad/sec) 

1 

9.870 

2 

39.478 

3 

88.825 

4 

157.912 

5 

246.730 

6 

355.301 

7 

483.605 

8 

631.647 

9 

799.428 

10 

986.950 

Table  3:  Mode  frequencies  for  pinned-pinned  beam 

not  contribute  to  G’(s).  From  (265)-(266)  it  can  be  shown  that  the  poles  in  Hbc  occur  at 

Sk  =  -$k2L2n2  ±  iyjl  -  ^k2L2-n2 

for  k  =  0,  ±1,  ±2, ....  For  these  computations  we  have  used  f  =  0.005  and  L  —  1.  The 
first  ten  modal  frequencies  are  shown  in  Table  2.  The  odd  numbered  (the  anti-symmetric) 
modes  only  appear  in  the  spectral  factor  as  can  be  seen  in  the  figures. 

The  distributed  gains  are  computed  according  to  (22)  by  numerical  quadrature  (trape¬ 
zoidal  rule).  The  resulting  gains  are  displayed  in  Figure  5.3.3.  The  gains  are  discontinuous 
at  the  point  of  observation  due  to  the  discontinuity  in  the  gradient  of  the  Green’s  function. 
A  major  issue  in  this  computation  is  the  numerical  precision  obtainable  by  integrating 
over  a  finite  bandwidth.  The  bandwidth  used  for  this  calculation  is  wo  =  1000 — the  same 
as  used  for  spectral  factorization.  Examination  of  the  magnitude  of  the  resulting  spec¬ 
tral  factor  indicates  that  with  the  material  damping  model  used  that  this  bandwidth  is 
reasonable  in  that  the  10‘h  mode  is  just  descernable  (in  double  precision)  in  the  spectral 
factor.  Examination  of  the  phase  response  indicates  however  that  significant  phase  excur¬ 
sions  are  still  obtained  at  high  frequencies  near  modal  frequencies.  This  phenomenon  can 
be  explained  in  terms  of  the  characteristic  interlacing  pattern  of  poles  and  zeros  for  this 
class  of  transfer  functions  [24,30].  Because  the  observation  point  z  =  L/2  is  not  ‘colocated’ 
with  the  control  point  z  =  L  the  transfer  function  will  have  nonminimum  phase  behavior. 
This  seems  to  lead  to  poor  convergence  of  the  distributed  gain  computations.  Figure  5.3.3 
displays  the  convergence  of  two  selected  points  of  the  distributed  gain 

1  0 

Kdi,{w )  =  —  /  (F*(tw)]-1G*(tw)CGr(iw;z,u;)  dw 

27T  J—u) o 

with  bandwidth  u>0.  Clearly  convergence  is  poor  as  the  bandwidth  of  the  integration  is 
increased.  The  effect  of  the  rapidly  varying  phase  is  apparently  significant. 

The  final  test  is  to  essentially  ‘re-compute’  the  achieved  spectral  factor  with  the  com¬ 
puted  Kn,.  This — in  effect — simulates  a  situation  where  the  distributed  state  feedback 
is  implemented  by  an  array  of  (in  this  case  100)  sensors  uniformly  distributed  across  the 
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Figure  11:  ideal  F{ iw)  -  Nyquiat  plot 


beam  measuring  both  components  of  the  state.  The  feedback  control  law  is  computed 
by  numerical  quadrature.  It  is  convenient  to  compare  the  spectral  factor  F(tu’)  com¬ 
puted  above  and  displayed  in  Figure  5.3.3  with  the  ‘achieved  return  difference’  F.(»u;)  = 
1  +  Jo  K4,,(w)Hbc(1u,  w)dw  displayed  in  Figure  5.3.3.  Figure  5.3.3  shows  the  character¬ 
istic  behavior  of  the  ideal  spectral  factor  with  respect  to  the  critical  point  at  a  =  0.  The 
Nyquist  contour  avoids  the  unit  disk  centered  at  the  origin.  In  Figure  5.3.3  the  achieved 
Nyquist  contour  indicates  a  phase  error  at  low  frequencies  and  for  higher  frequencies  it  ap¬ 
proach*  the  critical  point  rather  closely.  The  resulting  control  loop  will  be  only  marginally 
stable,  will  sustain  a  highly  underdamped  oscillation  at  some  high  fequency  and  will  be 
extremely  sensitive  to  modeling  errors  and  component  changes. 

At  this  point  the  above  negative  result  seems  to  indicate  a  fundamental  numerical 
sensitivity  associated  with  the  choice  of  ‘point’  control  as  an  optimzation  objective.  It 
seems  clear  that  by  projecting  the  cost  onto  a  finite  dimensional  ROM  subspace  will  lead 
to  a  much  more  tractable  problem  since  the  integration  can  now  be  performed  over  a 
known  finite  bandwidth  thus  assuring  convergence.  Of  Bourse  as  discussed  in  Section  2 
this  is  equivalent  to  solving  a  finite-dimensional  control  problem  for  the  ROM  model  so 
that  standard  methods  could  be  applied. 


£ 


6  Conclusions  and  Directions 


We  have  considered  modeling  of  flexible  structures  for  the  purpose  of  computing  con¬ 
trol  laws  for  distributed  parameter  systems.  We  have  considered  both  classical  modeling 
approachs  based  on  frequency  domain  analysis  for  well  defined  elastic  components  and 
homogenization  for  obtaining  effective  continuum  models  for  periodic  structures. 

A  complete  method  for  computation  of  distributed  control  laws  based  on  a  standard 
Wiener-Hopf  problem  has  been  demonstrated.  The  method  can,  in  principal,  provide 
optimal,  distributed  state  feedback  control  laws  for  systems  which  are  open  loop  stable.  In 
practice,  numerical  and  other  computational  considerations  serve  to  limit  the  effectiveness 
of  the  method  to  problems  with  well  defined  bandwidth  and  damping  which  provides 
a  uniform  exponential  rate  of  decay  for  essentially  all  modes.  Although  methods  were 
developed  which  can  reliably  determine  appropriate  irrational  transfer  function  models  for 
composite  structures  it  was  seen  that  computation  of  the  system  frequency  response  from 
such  transfer  functions  can  quickly  suffer  from  numerical  instability  outside  of  a  rather 
limited  bandwidth. 

It  was  shown  that  the  proposed  method  embraces  the  standard  notion  of  distributed 
control  by  reduced-order  approximation  by  definition  of  ROM  control  objective.  For  this 
case  the  major  advantage  of  the  method  appears  to  be  the  numerical  approach  which 
avoids  the  solution  (at  least  explicitly)  of  large  dimensional  Riccati  equation.  Of  course, 
as  discussed  by  Davis  [16],  it  is  often  desired  to  penalyze  only  those  modes  which  are 
physically  significant  for  the  design  at  hand.  In  our  experience  is  ultimately  required  to 
make  such  a  restriction  for  the  purposes  of  the  numerical  computations.  The  method  of 
computing  distributed  gains  is  itself  rather  sensitive  and  apparently  is  limited  in  practice 
to  computing  control  laws  for  only  a  finite  number  of  modes. 

The  frequency  response  method  is  particularly  appropriate  for  control  system  design 
and  system  performance  can  be  ascertained  in  advance  without  recourse  to  simulation.  By 
comparison  of  the  ideal  Nyquist  contour  for  the  optimal  problem  with  various  potential 
implimentations  of  distributed  state  sensing  resulting  stability  and  robustness  properties 
can  be  seen  graphically.  A  major  difficulty  with  practical  design  of  distributed  state  feed¬ 
back  is  that  an  array  of  sensors  spatially  separated  from  the  control  points  will  ultimately 
be  desired.  Current  advances  in  technology  offer  several  opportunities  for  relatively  low 
cost  alternatives  for  such  systems.  However,  whenever  sensors  are  spatially  separated  from 
actuators  for  such  systems  the  resulting  transfer  functions  are  nonminimum  phase.  Such 
an  effect  is  definitely  configuration  dependent  and  results  from  actuator  and  sensor  posi¬ 
tioning  but  has  a  dominant  effect  on  the  ability  to  provide  robust  and  high  performance 
control  laws. 

It  is  significant  to  note  the  success  of  the  spectral  factorization  computation  for  the 
problems  considered.  The  algorithm  is  quite  general  and  includes  the  case  of  multiple 
controls  where  the  required  factorization  is  for  a  matrix  function.  Application  of  this 
approach  to  the  computation  of  output  feedback  control  laws  for  distributed  systems  is 
considered  in  detail  in  [24],  Unlike  the  case  of  finite  dimensional  systems  the  output 
feedback  design  methods  we  have  tested  do  not  involve  explicit  distributed  state  estimation 
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from  sensor  measurements.  This  is  a  difficult  problem  which  definitely  requires  ROM 
approximation  for  on-line  computation.  Instead  we  make  use  of  alternate  state  realizations 
for  a  given  transfer  function  which  in  many  cases  consist  of  delay-differential  equations 
rather  than  partial  differential  equations.  These  equations  are  quite  easy  to  include  in  real¬ 
time  control  software.  Alternately,  if  such  realizations  cannot  be  found  we  can  often  realize 
the  control  by  general  purpose,  real-time  convolution.  Such  an  approach  is  particularly 
significant  given  the  increasing  availablity  of  extremely  high  speed  special  purpose  LSI 
chips  for  signal  processing.  Examples  of  these  include  Digital  Signal  Processing  (DSP) 
chips  and  Fast  Fourier  Transform  (FFT)  chips. 
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A  Algorithms  and  Code  for  Symbolic  Manipulation 

A  central  goal  of  this  project  was  to  investigate  the  use  of  a  symbolic  manipulation  (com¬ 
puter  algebra)  language  to  support  the  required  model  building  computations.  We  used  the 
language  SMP  throughout  this  project.  The  goal  was  to  automate — as  much  as  possible — 
computation  of  the  required  frequency  response  models  and  to  setup  the  equations  in  a 
form  suitable  numerical  evaluation.  Finally,  the  language  should  support  the  automatic 
generation  of  Fortran  code  which  can  then  be  linked  with  the  required  procedures  for 
spectral  factorization  and  optimal  gain  computation.  The  main  feature  of  the  algebraic 
approach  is  to  evaluate  the  expressions  symbolically  and  therefore  carry  through  certain 
model  parameters.  This  feature  is  quite  useful  for  generic  model  building  but  there  is  a 
penalty  in  that  complexity  of  the  resulting  expressions  can  limit  their  utility  while  limiting 
the  computational  speed  and  effectiveness  of  the  computer  algebra  system. 

The  first  step  was  to  build  a  system  shell  in  SMP  with  menu-driven,  user  interface 
which  can  aid  the  novice  user.  The  more  experienced  user  could  then  access  the  required 
SMP  procedures  directly.  As  with  any  computer  algebra  system  proficiency  comes  only 
with  practice  and  many  seemingly  straightforward  computations  can  lead  to  undesirably 
large  expressions.  Simplification  is  often  tedious  requiring  not  only  skill  with  the  special 
algebraic  relations  required  but  also  the  language  syntax  and  operations  for  symbolic  iden¬ 
tification  of  parts  of  expressions.  The  system  shell  we  implemented  included  the  equations 
for  computation  of  the  boundary  transfer  function  and  Green’s  function  for  both  the  hy¬ 
perbolic  and  parabolic  forms  developed  in  Section  3.  We  also  included  the  equations  for 
computation  of  hybrid  coordinate  models. 

An  essential  computation  for  these  equations  is  evaluation  of  a  matrix  exponential.  We 
found  it  convenient  to  use  an  algebraic  relation  based  on  the  Laplace  transform  and  the 
Cayley-Hamilton  theorem  for  matrices.  The  relation  is  commonly  credited  to  Leverrier- 
Souriau-Faddeeva- Frame  [85,  pp  658].  The  formulae  are  developed  as  follows.  Given  an 
n  x  n  matrix  A  the  characteristic  equation  is 

xa(s)  =  det[s/  -  A] 

=  oq  +  oi«  +  . . .  +  on_ isn  1  +  s". 

Then  the  Cayley-Hamilton  theorem  says  that 


*a{A)  =  0. 


The  matrix  exponential 


can  then  be  written  as  a  finite  sum; 


k=0 
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The  coefficients  can  be  computed  as 

iMO  f  ~—fT e‘ds 

2jrt  Jr  *,*(«) 

for  T  a  smooth,  closed  contour  enclosing  the  spectrum  of  A  and  where  the  polynomials 
jrM(s)  are  generated  recursively  according  to 

jrl,-1J(s)  =  «Jr(i)(s)  +  a/_iJr(n)(s) 


with 

*(n)(s)  =  i. 

The  contour  integral  is  the  usual  Laplace  transform  inverse  which  can  be  evaluated  by 
residues.  Most  available  computer  algebra  systems  (e.g.  MACSYMA,  SMP,  REDUCE) 
can  perform  the  calculation  by  explicit  partial  fraction  expansion  of  the  rational  integrand 
in  combination  with  certain  tables  for  standard  inverses.  This  approach  was  programmed 
in  SMP  and  tested.  It  proved  considerably  faster  than  the  direct  approach;  i.e.,  compute 
eAt  by  inverse  Laplace  transform  of  the  n  x  n  matrix  [sJ  —  A]-1. 

SMP  offers  the  facility  for  ‘tuning’  automatic  simplification  of  expressions  by  adding 
‘rules’  which  are  automatically  applied  at  each  computation.  We  investigated  this  feature 
for  automatic  simplication  of  the  transcendental  terms  generated  for  transfer  function  com¬ 
putation.  Unfortunately  it  proved  to  be  difficult  to  determine  a  set  of  algebraic  rules  which 
would  always  lead  to  “desirable”  reduced  forms.  An  alternate  approach  we  found  useful 
was  to  build  special  simplification  procedures  in  SMP  for  various  steps  in  the  computation. 
These  procedures  could  be  applied  manually  to  resulting  expressions.  They  could  also  be 
applied  with  the  *MAP[]’  operation  of  SMP  which  recursively  applies  a  set  of  procedures 
until  the  expression  does  not  change.  This  can  lead  to  indefinite  recursions  for  certain 
expressions.  In  conclusion,  we  found  that  SMP  offered  several  features  which  were  useful 
for  expression  simplification  but  no  automatic  procedure  could  be  developed  which  would 
always  lead  to  a  ‘nice  enough’  expression  that  the  knowledgable  human  could  not  ask  for 
more.  In  the  final  analysis  computer  algebra  is  a  tool  for  experts— not  (in  the  current 
jargon)  an  expert  system! 

One  aspect  of  computer  algebra  that  we  found  most  useful  on  this  project  however  is  the 
use  of  standard  expansions  for  generating  alternate  expressions  for  numerically  evaluating 
the  required  irrational  transfer  functions.  This  was  an  essential  task  and  the  understanding 
of  the  numerical  stability  properties  of  transcendental  terms  absorbed  a  major  portion  of 
this  project.  In  order  to  make  a  computer  system  ultimately  useful  to  engineers  it  will 
be  necessary  to  include  specific  provisions  for  numerical  testing  and  analysis.  Coupling 
this  with  a  feature  for  automatic  Fortran  or  C  code  generation  should  provide  a  useful 
tool  for  advanced  model  development  and  analysis  of  mechanical  systems.  The  feature  for 
automatic  Fortran  code  generation  was  not  fully  debugged  in  the  version  of  SMP  that  we 
used.  We  were  able  to  compensate  for  most  of  the  bugs  with  a  programmable  editor  (Gnu 
Emacs). 

The  following  code  is  the  SMP  system  shell  we  developed  for  this  project.  It  includes  a 
menu-driven,  user  interface  containing  all  the  required  equations  for  computing  the  transfer 
function  models  for  hyperbolic,  parabolic,  and  hybrid  configurations. 


Menu  Driven  Programs 


Main  Menu 

Distributed  Parameter  System  (DPS) 
Lumped  Parameter  System  (LPS) 
Hybrid  System 
Parabolic  System 
Incidence  Relations 
Cost  Definition 
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Matrix  Exponential 
Augment  Matrices 
Join  Matrices 
Block  Diagonalize  Matrix 
Extract  Matrix 
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Miscellaneous  Routines 


Power  Series  Expansion 
Reduce  subexpressions 
Get  external  files 
Define  properties  to  symbols 
Generate  FORTRAN  CODE 
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REAL  FUNCTION  RFUNC1 (ABC ,LC1 ,LC2 ,LC3,KBYT, INC , npt*) 
•••  INSERT  COMMON  BLOCK  FOR  STACK  HERE  ••• 

INCLUDE  ’SPEC  STACK. inc' 

INTEGER  KABC,RORO,NOIM,INCl 
DATA  KABC,KORO,NDIM,INCl/0, 0,0,0/ 
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COPIES  REAL/ IMAGINARY  PARTS  OP  A  VECTOR  INTO  COMPLEX  ARRAY 
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